commit 7f9396087b7378a3afce69c99e3f64465b4af0e7
parent 3c3aa3031dee8ce14cd68d78c0092d6900dfb669
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Tue, 5 May 2020 22:27:26 +0200
Use consistent formatting style
Diffstat:
| M | 1d_fd_simple_shear.c | | | 25 | +++++++++++++------------ |
| M | arrays.c | | | 153 | ++++++++++++++++++++++++++++++++++++++++++++----------------------------------- |
| M | fluid.c | | | 203 | +++++++++++++++++++++++++++++++++++++++++-------------------------------------- |
| M | max_depth_simple_shear.c | | | 96 | ++++++++++++++++++++++++++++++++++++++++--------------------------------------- |
| M | shear_flux.c | | | 4 | ++-- |
| M | simulation.c | | | 389 | +++++++++++++++++++++++++++++++++++++++---------------------------------------- |
6 files changed, 447 insertions(+), 423 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -14,7 +14,7 @@
#define MAX_ITER_1D_FD_SIMPLE_SHEAR 100000
/* uncomment to print time spent per time step to stdout */
-/*#define BENCHMARK_PERFORMANCE*/
+/* #define BENCHMARK_PERFORMANCE */
char *argv0;
@@ -63,15 +63,17 @@ usage(void)
}
int
-main(int argc, char* argv[])
+main(int argc, char *argv[])
{
int i, normalize;
unsigned long iter;
double new_phi, new_k, filetimeclock;
struct simulation sim;
+
#ifdef BENCHMARK_PERFORMANCE
clock_t t_begin, t_end;
double t_elapsed;
+
#endif
#ifdef __OpenBSD__
@@ -180,7 +182,8 @@ main(int argc, char* argv[])
argv[1]);
return 1;
}
- argc--; argv++;
+ argc--;
+ argv++;
break;
case 's':
sim.v_x_fix = atof(EARGF(usage()));
@@ -201,7 +204,7 @@ main(int argc, char* argv[])
sim.v_x_bot = atof(EARGF(usage()));
break;
case 'v':
- printf("%s-"VERSION"\n", argv0);
+ printf("%s-" VERSION "\n", argv0);
return 0;
case 'Y':
sim.phi_max = atof(EARGF(usage()));
@@ -209,7 +212,7 @@ main(int argc, char* argv[])
case 'y':
sim.phi_min = atof(EARGF(usage()));
break;
- default:
+default:
usage();
} ARGEND;
@@ -219,15 +222,15 @@ main(int argc, char* argv[])
usage();
if (sim.nz < 1)
- sim.nz = (int)ceil(sim.L_z/sim.d);
+ sim.nz = (int) ceil(sim.L_z / sim.d);
prepare_arrays(&sim);
if (!isnan(new_phi))
- for (i=0; i<sim.nz; ++i)
+ for (i = 0; i < sim.nz; ++i)
sim.phi[i] = new_phi;
if (!isnan(new_k))
- for (i=0; i<sim.nz; ++i)
+ for (i = 0; i < sim.nz; ++i)
sim.k[i] = new_k;
lithostatic_pressure_distribution(&sim);
@@ -257,20 +260,18 @@ main(int argc, char* argv[])
free_arrays(&sim);
exit(10);
}
-
#ifdef BENCHMARK_PERFORMANCE
t_end = clock();
- t_elapsed = (double)(t_end - t_begin)/CLOCKS_PER_SEC;
+ t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC;
printf("time spent per time step = %g s\n", t_elapsed);
#endif
if ((filetimeclock >= sim.file_dt || iter == 1) &&
strncmp(sim.name, DEFAULT_SIMULATION_NAME,
- sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
+ sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
write_output_file(&sim, normalize);
filetimeclock = 0.0;
}
-
sim.t += sim.dt;
filetimeclock += sim.dt;
iter++;
diff --git a/arrays.c b/arrays.c
@@ -5,17 +5,17 @@
#define DEG2RAD(x) (x*M_PI/180.0)
void
-check_magnitude(const char* func_name, int limit, int value)
+check_magnitude(const char *func_name, int limit, int value)
{
if (value < limit) {
fprintf(stderr, "error: %s: input size %d is less than %d\n",
- func_name, value, limit);
+ func_name, value, limit);
exit(1);
}
}
-/* Translate a i,j,k index in grid with dimensions nx, ny, nz into a linear
- * index */
+/* Translate a i,j,k index in grid with dimensions nx, ny, nz into a
+ * linear index */
unsigned int
idx3(const unsigned int i,
const unsigned int j,
@@ -23,11 +23,11 @@ idx3(const unsigned int i,
const unsigned int nx,
const unsigned int ny)
{
- return i + nx*j + nx*ny*k;
+ return i + nx * j + nx * ny * k;
}
-/* Translate a i,j,k index in grid with dimensions nx, ny, nz and a padding of
- * single ghost nodes into a linear index */
+/* Translate a i,j,k index in grid with dimensions nx, ny, nz and a
+ * padding of single ghost nodes into a linear index */
unsigned int
idx3g(const unsigned int i,
const unsigned int j,
@@ -35,33 +35,36 @@ idx3g(const unsigned int i,
const unsigned int nx,
const unsigned int ny)
{
- return i+1 + (nx+2)*(j+1) + (nx+2)*(ny+2)*(k+1);
+ return i + 1 + (nx + 2) * (j + 1) + (nx + 2) * (ny + 2) * (k + 1);
}
-/* Translate a i,j index in grid with dimensions nx, ny into a linear index */
+/* Translate a i,j index in grid with dimensions nx, ny into a linear
+ * index */
unsigned int
idx2(const unsigned int i, const unsigned int j, const unsigned int nx)
{
- return i + nx*j;
+ return i + nx * j;
}
-/* Translate a i,j index in grid with dimensions nx, ny and a padding of single
- * ghost nodes into a linear index */
+/* Translate a i,j index in grid with dimensions nx, ny and a padding of
+ * single ghost nodes into a linear index */
unsigned int
idx2g(const unsigned int i, const unsigned int j, const unsigned int nx)
{
- return i+1 + (nx+2)*(j+1);
+ return i + 1 + (nx + 2) * (j + 1);
}
-/* Translate a i index in grid with a padding of single into a linear index */
+/* Translate a i index in grid with a padding of single into a linear
+ * index */
unsigned int
idx1g(const unsigned int i)
{
- return i+1;
+ return i + 1;
}
-/* Return an array of `n` linearly spaced values in the range [lower; upper] */
-double*
+/* Return an array of `n` linearly spaced values in the range [lower;
+ * upper] */
+double *
linspace(const double lower, const double upper, const int n)
{
int i;
@@ -69,205 +72,221 @@ linspace(const double lower, const double upper, const int n)
double dx;
check_magnitude(__func__, 1, n);
- x = malloc(n*sizeof(double));
- dx = (upper - lower)/(double)(n-1);
- for (i=0; i<n; ++i)
- x[i] = lower + dx*i;
+ x = malloc(n * sizeof(double));
+ dx = (upper - lower) / (double) (n - 1);
+ for (i = 0; i < n; ++i)
+ x[i] = lower + dx * i;
+
return x;
}
/* Return an array of `n-1` values with the intervals between `x` values */
-double*
-spacing(const double* x, const int n)
+double *
+spacing(const double *x, const int n)
{
int i;
double *dx;
check_magnitude(__func__, 2, n);
- dx = malloc((n-1)*sizeof(double));
- for (i=0; i<n-1; ++i)
- dx[i] = x[i+1] - x[i];
+ dx = malloc((n - 1) * sizeof(double));
+ for (i = 0; i < n - 1; ++i)
+ dx[i] = x[i + 1] - x[i];
+
return dx;
}
/* Return an array of `n` values with the value 0.0 */
-double*
+double *
zeros(const int n)
{
int i;
double *x;
check_magnitude(__func__, 1, n);
- x = malloc(n*sizeof(double));
- for (i=0; i<n; ++i)
+ x = malloc(n * sizeof(double));
+ for (i = 0; i < n; ++i)
x[i] = 0.0;
+
return x;
}
/* Return an array of `n` values with the value 1.0 */
-double*
+double *
ones(const int n)
{
int i;
double *x;
check_magnitude(__func__, 1, n);
- x = malloc(n*sizeof(double));
- for (i=0; i<n; ++i)
+ x = malloc(n * sizeof(double));
+ for (i = 0; i < n; ++i)
x[i] = 1.0;
+
return x;
}
/* Return an array of `n` values with a specified value */
-double*
+double *
initval(const double value, const int n)
{
int i;
double *x;
check_magnitude(__func__, 1, n);
- x = malloc(n*sizeof(double));
- for (i=0; i<n; ++i)
+ x = malloc(n * sizeof(double));
+ for (i = 0; i < n; ++i)
x[i] = value;
+
return x;
}
/* Return an array of `n` uninitialized values */
-double*
+double *
empty(const int n)
{
check_magnitude(__func__, 1, n);
- return malloc(n*sizeof(double));
+ return malloc(n * sizeof(double));
}
/* Return largest value in array `a` with size `n` */
double
-max(const double* a, const int n)
+max(const double *a, const int n)
{
int i;
double maxval;
check_magnitude(__func__, 1, n);
maxval = -INFINITY;
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
if (a[i] > maxval)
maxval = a[i];
+
return maxval;
}
/* Return smallest value in array `a` with size `n` */
double
-min(const double* a, const int n)
+min(const double *a, const int n)
{
int i;
double minval;
check_magnitude(__func__, 1, n);
minval = +INFINITY;
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
if (a[i] < minval)
minval = a[i];
+
return minval;
}
void
-print_array(const double* a, const int n)
+print_array(const double *a, const int n)
{
int i;
+
check_magnitude(__func__, 1, n);
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
printf("%.17g\n", a[i]);
}
void
-print_arrays(const double* a, const double* b, const int n)
+print_arrays(const double *a, const double *b, const int n)
{
int i;
+
check_magnitude(__func__, 1, n);
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
printf("%.17g\t%.17g\n", a[i], b[i]);
}
void
-print_arrays_2nd_normalized(const double* a, const double* b, const int n)
+print_arrays_2nd_normalized(const double *a, const double *b, const int n)
{
int i;
double max_b;
+
check_magnitude(__func__, 1, n);
max_b = max(b, n);
- for (i=0; i<n; ++i)
- printf("%.17g\t%.17g\n", a[i], b[i]/max_b);
+ for (i = 0; i < n; ++i)
+ printf("%.17g\t%.17g\n", a[i], b[i] / max_b);
}
void
-print_three_arrays(const double* a,
- const double* b,
- const double* c,
+print_three_arrays(const double *a,
+ const double *b,
+ const double *c,
const int n)
{
int i;
+
check_magnitude(__func__, 1, n);
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
}
void
-fprint_arrays(FILE* fp, const double* a, const double* b, const int n)
+fprint_arrays(FILE *fp, const double *a, const double *b, const int n)
{
int i;
+
check_magnitude(__func__, 1, n);
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]);
}
void
-fprint_three_arrays(FILE* fp,
- const double* a,
- const double* b,
- const double* c,
+fprint_three_arrays(FILE *fp,
+ const double *a,
+ const double *b,
+ const double *c,
const int n)
{
int i;
+
check_magnitude(__func__, 1, n);
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]);
}
void
-copy_values(const double* in, double* out, const int n)
+copy_values(const double *in, double *out, const int n)
{
int i;
+
check_magnitude(__func__, 1, n);
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
out[i] = in[i];
}
-double*
-copy(const double* in, const int n)
+double *
+copy(const double *in, const int n)
{
double *out;
+
check_magnitude(__func__, 1, n);
out = empty(n);
copy_values(in, out, n);
return out;
}
-double*
-normalize(const double* in, const int n)
+double *
+normalize(const double *in, const int n)
{
int i;
double max_val;
double *out;
check_magnitude(__func__, 1, n);
- out = malloc(n*sizeof(double));
+ out = malloc(n * sizeof(double));
copy_values(in, out, n);
max_val = max(out, n);
if (max_val == 0.0)
max_val = 1.0;
- for (i=0; i<n; ++i)
+ for (i = 0; i < n; ++i)
out[i] /= max_val;
+
return out;
}
diff --git a/fluid.c b/fluid.c
@@ -9,16 +9,17 @@ void
hydrostatic_fluid_pressure_distribution(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->p_f_ghost[i+1] = sim->p_f_top
- + sim->phi[i]*sim->rho_f*sim->G
- *(sim->L_z - sim->z[i]);
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->p_f_ghost[i + 1] = sim->p_f_top
+ + sim->phi[i] * sim->rho_f * sim->G
+ * (sim->L_z - sim->z[i]);
}
-/* Determines the largest time step for the current simulation state. Note
- * that the time step should be recalculated if cell sizes or diffusivities
- * (i.e., permeabilities, porosities, viscosities, or compressibilities)
- * change. The safety factor should be in ]0;1] */
+/* Determines the largest time step for the current simulation state. Note
+ * that the time step should be recalculated if cell sizes or
+ * diffusivities (i.e., permeabilities, porosities, viscosities, or
+ * compressibilities) change. The safety factor should be in ]0;1] */
int
set_largest_fluid_timestep(struct simulation *sim, const double safety)
{
@@ -28,26 +29,29 @@ set_largest_fluid_timestep(struct simulation *sim, const double safety)
dx = spacing(sim->z, sim->nz);
dx_min = INFINITY;
- for (i=0; i<sim->nz-1; ++i) {
+ for (i = 0; i < sim->nz - 1; ++i) {
if (dx[i] < 0.0) {
fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n",
dx[i], i);
free(dx);
return 1;
}
- if (dx[i] < dx_min) dx_min = dx[i];
+ if (dx[i] < dx_min)
+ dx_min = dx[i];
}
free(dx);
diff_max = -INFINITY;
- for (i=0; i<sim->nz; ++i) {
- diff = sim->k[i]/(sim->phi[i]*sim->beta_f*sim->mu_f);
- if (diff > diff_max) diff_max = diff;
+ for (i = 0; i < sim->nz; ++i) {
+ diff = sim->k[i] / (sim->phi[i] * sim->beta_f * sim->mu_f);
+ if (diff > diff_max)
+ diff_max = diff;
}
- sim->dt = safety*0.5*dx_min*dx_min/diff_max;
- if (sim->file_dt*0.5 < sim->dt)
+ sim->dt = safety * 0.5 * dx_min * dx_min / diff_max;
+ if (sim->file_dt * 0.5 < sim->dt)
sim->dt = sim->file_dt;
+
return 0;
}
@@ -57,7 +61,7 @@ sine_wave(const double time,
const double freq,
const double phase)
{
- return ampl*sin(2.0*PI*freq*time + phase);
+ return ampl * sin(2.0 * PI * freq * time + phase);
}
static double
@@ -66,10 +70,10 @@ triangular_pulse(const double time,
const double freq,
const double peak_time)
{
- if (peak_time - 1.0/(2.0*freq) < time && time <= peak_time)
- return peak_ampl*2.0*freq*(time - peak_time) + peak_ampl;
- else if (peak_time < time && time < peak_time + 1.0/(2.0*freq))
- return peak_ampl*2.0*freq*(peak_time - time) + peak_ampl;
+ if (peak_time - 1.0 / (2.0 * freq) < time && time <= peak_time)
+ return peak_ampl * 2.0 * freq * (time - peak_time) + peak_ampl;
+ else if (peak_time < time && time < peak_time + 1.0 / (2.0 * freq))
+ return peak_ampl * 2.0 * freq * (peak_time - time) + peak_ampl;
else
return 0.0;
}
@@ -80,8 +84,8 @@ square_pulse(const double time,
const double freq,
const double peak_time)
{
- if (peak_time - 1.0/(2.0*freq) < time &&
- time < peak_time + 1.0/(2.0*freq))
+ if (peak_time - 1.0 / (2.0 * freq) < time &&
+ time < peak_time + 1.0 / (2.0 * freq))
return peak_ampl;
else
return 0.0;
@@ -91,52 +95,58 @@ static void
set_fluid_bcs(struct simulation *sim, const double p_f_top)
{
set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top);
- sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */
+ sim->p_f_ghost[sim->nz] = p_f_top; /* Include top node in BC */
set_bc_neumann(sim->p_f_ghost,
sim->nz,
-1,
- sim->phi[0]*sim->rho_f*sim->G,
+ sim->phi[0] * sim->rho_f * sim->G,
sim->dz);
}
static double
darcy_pressure_change_1d(const int i,
- const int nz,
- const double* p_f_ghost_in,
- const double* phi,
- const double* phi_dot,
- const double* k,
- const double dz,
- const double beta_f,
- const double mu_f)
+ const int nz,
+ const double *p_f_ghost_in,
+ const double *phi,
+ const double *phi_dot,
+ const double *k,
+ const double dz,
+ const double beta_f,
+ const double mu_f)
{
double k_ = k[i], div_k_grad_p, k_zn, k_zp;
- if (i==0) k_zn = k_; else k_zn = k[i-1];
- if (i==nz-1) k_zp = k_; else k_zp = k[i+1];
+ if (i == 0)
+ k_zn = k_;
+ else
+ k_zn = k[i - 1];
+ if (i == nz - 1)
+ k_zp = k_;
+ else
+ k_zp = k[i + 1];
#ifdef DEBUG
printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n",
- i, i+1,
- p_f_ghost_in[i], p_f_ghost_in[i+1], p_f_ghost_in[i+2],
- k_zn, k_, k_zp);
+ i, i + 1,
+ p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
+ k_zn, k_, k_zp);
#endif
- div_k_grad_p = (2.0*k_zp*k_/(k_zp + k_)
- * (p_f_ghost_in[i+2]
- - p_f_ghost_in[i+1])/dz
- - 2.0*k_zn*k_/(k_zn + k_)
- * (p_f_ghost_in[i+1]
- - p_f_ghost_in[i])/dz
- )/dz;
+ div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
+ * (p_f_ghost_in[i + 2]
+ - p_f_ghost_in[i + 1]) / dz
+ - 2.0 * k_zn * k_ / (k_zn + k_)
+ * (p_f_ghost_in[i + 1]
+ - p_f_ghost_in[i]) / dz
+ ) / dz;
#ifdef DEBUG
printf("darcy_pressure_change [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
- i, phi[i], div_k_grad_p, phi_dot[i]);
+ i, phi[i], div_k_grad_p, phi_dot[i]);
#endif
/* TODO: add advective term */
- return 1.0/(beta_f*phi[i]*mu_f)*div_k_grad_p
- - 1.0/(beta_f*(1.0 - phi[i]))*phi_dot[i];
+ return 1.0 / (beta_f * phi[i] * mu_f) * div_k_grad_p
+ - 1.0 / (beta_f * (1.0 - phi[i])) * phi_dot[i];
}
int
@@ -148,39 +158,36 @@ darcy_solver_1d(struct simulation *sim,
double epsilon, theta, p_f_top, r_norm_max = NAN;
double *dp_f_dt_expl;
double *p_f_ghost_old, *dp_f_dt_impl, *p_f_ghost_new, *r_norm;
-
- /* choose integration method, parameter in [0.0; 1.0]
- * epsilon = 0.0: explicit
- * epsilon = 0.5: Crank-Nicolson
- * epsilon = 1.0: implicit */
+
+ /* choose integration method, parameter in [0.0; 1.0] epsilon =
+ * 0.0: explicit epsilon = 0.5: Crank-Nicolson epsilon = 1.0:
+ * implicit */
epsilon = 0.5;
- /* choose relaxation factor, parameter in ]0.0; 1.0]
- * theta in ]0.0; 1.0]: underrelaxation
- * theta = 1.0: Gauss-Seidel
- * theta > 1.0: overrelaxation */
+ /* choose relaxation factor, parameter in ]0.0; 1.0] theta in
+ * ]0.0; 1.0]: underrelaxation theta = 1.0: Gauss-Seidel theta >
+ * 1.0: overrelaxation */
/* TODO: values other than 1.0 do not work! */
theta = 1.0;
if (isnan(sim->p_f_mod_pulse_time))
- p_f_top = sim->p_f_top
- + sine_wave(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_phase);
+ p_f_top = sim->p_f_top
+ + sine_wave(sim->t,
+ sim->p_f_mod_ampl,
+ sim->p_f_mod_freq,
+ sim->p_f_mod_phase);
+ else if (sim->p_f_mod_pulse_shape == 1)
+ p_f_top = sim->p_f_top
+ + square_pulse(sim->t,
+ sim->p_f_mod_ampl,
+ sim->p_f_mod_freq,
+ sim->p_f_mod_pulse_time);
else
- if (sim->p_f_mod_pulse_shape == 1)
- p_f_top = sim->p_f_top
- + square_pulse(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_pulse_time);
- else
- p_f_top = sim->p_f_top
- + triangular_pulse(sim->t,
- sim->p_f_mod_ampl,
- sim->p_f_mod_freq,
- sim->p_f_mod_pulse_time);
+ p_f_top = sim->p_f_top
+ + triangular_pulse(sim->t,
+ sim->p_f_mod_ampl,
+ sim->p_f_mod_freq,
+ sim->p_f_mod_pulse_time);
/* set fluid BCs (1 of 2) */
set_fluid_bcs(sim, p_f_top);
@@ -188,7 +195,7 @@ darcy_solver_1d(struct simulation *sim,
if (epsilon < 1.0) {
/* compute explicit solution to pressure change */
dp_f_dt_expl = zeros(sim->nz);
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
dp_f_dt_expl[i] = darcy_pressure_change_1d(i,
sim->nz,
sim->p_f_ghost,
@@ -199,25 +206,24 @@ darcy_solver_1d(struct simulation *sim,
sim->beta_f,
sim->mu_f);
}
-
if (epsilon > 0.0) {
/* compute implicit solution with Jacobian iterations */
dp_f_dt_impl = zeros(sim->nz);
- p_f_ghost_new = zeros(sim->nz+2);
+ p_f_ghost_new = zeros(sim->nz + 2);
r_norm = zeros(sim->nz);
- p_f_ghost_old = empty(sim->nz+2);
- copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz+2);
+ p_f_ghost_old = empty(sim->nz + 2);
+ copy_values(sim->p_f_ghost, p_f_ghost_old, sim->nz + 2);
- for (iter=0; iter<max_iter; ++iter) {
+ for (iter = 0; iter < max_iter; ++iter) {
/* set fluid BCs (2 of 2) */
set_fluid_bcs(sim, p_f_top);
#ifdef DEBUG
puts(".. p_f_ghost after BC:");
- print_array(sim->p_f_ghost, sim->nz+2);
+ print_array(sim->p_f_ghost, sim->nz + 2);
#endif
- for (i=0; i<sim->nz-1; ++i)
+ for (i = 0; i < sim->nz - 1; ++i)
dp_f_dt_impl[i] = darcy_pressure_change_1d(i,
sim->nz,
sim->p_f_ghost,
@@ -228,37 +234,37 @@ darcy_solver_1d(struct simulation *sim,
sim->beta_f,
sim->mu_f);
- for (i=0; i<sim->nz-1; ++i) {
+ for (i = 0; i < sim->nz - 1; ++i) {
#ifdef DEBUG
printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n",
- i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
+ i, dp_f_dt_expl[i], i, dp_f_dt_impl[i]);
#endif
/* TODO */
- p_f_ghost_new[i+1] = p_f_ghost_old[i+1]
- + epsilon*dp_f_dt_impl[i]*sim->dt;
+ p_f_ghost_new[i + 1] = p_f_ghost_old[i + 1]
+ + epsilon * dp_f_dt_impl[i] * sim->dt;
if (epsilon < 1.0)
- p_f_ghost_new[i+1] += (1.0 - epsilon)
- *dp_f_dt_expl[i]*sim->dt;
+ p_f_ghost_new[i + 1] += (1.0 - epsilon)
+ * dp_f_dt_expl[i] * sim->dt;
- p_f_ghost_new[i+1] = p_f_ghost_old[i+1]*(1.0 - theta)
- + p_f_ghost_new[i+1]*theta;
+ p_f_ghost_new[i + 1] = p_f_ghost_old[i + 1] * (1.0 - theta)
+ + p_f_ghost_new[i + 1] * theta;
- r_norm[i] = residual_normalized(p_f_ghost_new[i+1],
- sim->p_f_ghost[i+1]);
+ r_norm[i] = residual_normalized(p_f_ghost_new[i + 1],
+ sim->p_f_ghost[i + 1]);
}
r_norm_max = max(r_norm, sim->nz);
#ifdef DEBUG
puts(".. p_f_ghost_new:");
- print_array(p_f_ghost_new, sim->nz+2);
+ print_array(p_f_ghost_new, sim->nz + 2);
#endif
- copy_values(p_f_ghost_new, sim->p_f_ghost, sim->nz+2);
+ copy_values(p_f_ghost_new, sim->p_f_ghost, sim->nz + 2);
#ifdef DEBUG
puts(".. p_f_ghost after update:");
- print_array(sim->p_f_ghost, sim->nz+2);
+ print_array(sim->p_f_ghost, sim->nz + 2);
#endif
if (r_norm_max <= rel_tol) {
@@ -276,23 +282,24 @@ darcy_solver_1d(struct simulation *sim,
if (!solved) {
fprintf(stderr, "darcy_solver_1d: ");
fprintf(stderr, "Solution did not converge after %d iterations\n",
- iter);
+ iter);
fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
}
} else {
- for (i=0; i<sim->nz; ++i)
- sim->p_f_ghost[i+1] += dp_f_dt_expl[i]*sim->dt;
+ for (i = 0; i < sim->nz; ++i)
+ sim->p_f_ghost[i + 1] += dp_f_dt_expl[i] * sim->dt;
solved = 1;
#ifdef DEBUG
puts(".. dp_f_dt_expl:");
print_array(dp_f_dt_expl, sim->nz);
puts(".. p_f_ghost after explicit solution:");
- print_array(sim->p_f_ghost, sim->nz+2);
+ print_array(sim->p_f_ghost, sim->nz + 2);
#endif
}
set_fluid_bcs(sim, p_f_top);
if (epsilon < 1.0)
free(dp_f_dt_expl);
+
return solved - 1;
}
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -46,30 +46,30 @@ usage(void)
static double
skin_depth(const struct simulation *sim)
{
- return sqrt(sim->k[0]/
- (sim->phi[0]*sim->mu_f*sim->beta_f*M_PI*sim->p_f_mod_freq));
+ return sqrt(sim->k[0] /
+ (sim->phi[0] * sim->mu_f * sim->beta_f * M_PI * sim->p_f_mod_freq));
}
/* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */
static double
-eff_normal_stress_gradient(struct simulation *sim, double d_s, double z_)
+eff_normal_stress_gradient(struct simulation * sim, double d_s, double z_)
{
- if (z_/d_s > 10.0) {
+ if (z_ / d_s > 10.0) {
fprintf(stderr, "error: %s: unrealistic depth: %g m "
- "relative to skin depth %g m\n", __func__, z_, d_s);
+ "relative to skin depth %g m\n", __func__, z_, d_s);
free_arrays(sim);
exit(10);
}
-
- return sqrt(2.0)*sin((3.0*M_PI/2.0 - z_/d_s) + M_PI/4.0)
- + (sim->rho_s - sim->rho_f)*sim->G*d_s/sim->p_f_mod_ampl*exp(z_/d_s);
+ return sqrt(2.0) * sin((3.0 * M_PI / 2.0 - z_ / d_s) + M_PI / 4.0)
+ + (sim->rho_s - sim->rho_f) * sim->G * d_s
+ / sim->p_f_mod_ampl * exp(z_ / d_s);
}
/* use Brent's method for finding roots (Press et al., 2007) */
static double
zbrent(struct simulation *sim,
double d_s,
- double (*f)(struct simulation *sim, double, double),
+ double (*f) (struct simulation *sim, double, double),
double x1,
double x2,
double tol)
@@ -80,18 +80,18 @@ zbrent(struct simulation *sim,
a = x1;
b = x2;
c = x2;
- fa = (*f)(sim, d_s, a);
- fb = (*f)(sim, d_s, b);
+ fa = (*f) (sim, d_s, a);
+ fb = (*f) (sim, d_s, b);
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
fprintf(stderr,
- "error: %s: no root in range %g,%g m\n",
- __func__, x1, x2);
+ "error: %s: no root in range %g,%g m\n",
+ __func__, x1, x2);
free_arrays(sim);
exit(11);
}
fc = fb;
- for (iter=0; iter<MAX_ITER; iter++) {
+ for (iter = 0; iter < MAX_ITER; iter++) {
if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
c = a;
fc = fa;
@@ -106,29 +106,29 @@ zbrent(struct simulation *sim,
fb = fc;
fc = fa;
}
- tol1 = 2.0*EPS*fabs(b) + 0.5*tol;
- xm = 0.5*(c - b);
+ tol1 = 2.0 * EPS * fabs(b) + 0.5 * tol;
+ xm = 0.5 * (c - b);
if (fabs(xm) <= tol1 || fb == 0.0)
return b;
if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
- s = fb/fa;
+ s = fb / fa;
if (a == c) {
- p = 2.0*xm*s;
+ p = 2.0 * xm * s;
q = 1.0 - s;
} else {
- q = fa/fc;
- r = fb/fc;
- p = s*(2.0*xm*q*(q - r) - (b - a)*(r - 1.0));
- q = (q - 1.0)*(r - 1.0)*(s - 1.0);
+ q = fa / fc;
+ r = fb / fc;
+ p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
+ q = (q - 1.0) * (r - 1.0) * (s - 1.0);
}
if (p > 0.0)
q = -q;
p = fabs(p);
- min1 = 3.0*xm*q - fabs(tol1*q);
- min2 = fabs(e*q);
- if (2.0*p < (min1 < min2 ? min1 : min2)) {
+ min1 = 3.0 * xm * q - fabs(tol1 * q);
+ min2 = fabs(e * q);
+ if (2.0 * p < (min1 < min2 ? min1 : min2)) {
e = d;
- d = p/q;
+ d = p / q;
} else {
d = xm;
e = d;
@@ -143,26 +143,28 @@ zbrent(struct simulation *sim,
b += d;
else
b += ((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1));
- fb = (*f)(sim, d_s, b);
+ fb = (*f) (sim, d_s, b);
}
fprintf(stderr,
- "error: %s: exceeded maximum number of iterations",
- __func__);
+ "error: %s: exceeded maximum number of iterations",
+ __func__);
free_arrays(sim);
exit(12);
return NAN;
}
int
-main(int argc, char* argv[])
+main(int argc, char *argv[])
{
int i;
double new_phi, new_k;
double d_s, depth, depth_limit1, depth_limit2;
struct simulation sim;
+
#ifdef BENCHMARK_PERFORMANCE
clock_t t_begin, t_end;
double t_elapsed;
+
#endif
#ifdef __OpenBSD__
@@ -211,7 +213,7 @@ main(int argc, char* argv[])
sim.rho_s = atof(EARGF(usage()));
break;
case 'v':
- printf("%s-"VERSION"\n", argv0);
+ printf("%s-" VERSION "\n", argv0);
return 0;
break;
default:
@@ -225,10 +227,10 @@ main(int argc, char* argv[])
prepare_arrays(&sim);
if (!isnan(new_phi))
- for (i=0; i<sim.nz; ++i)
+ for (i = 0; i < sim.nz; ++i)
sim.phi[i] = new_phi;
if (!isnan(new_k))
- for (i=0; i<sim.nz; ++i)
+ for (i = 0; i < sim.nz; ++i)
sim.k[i] = new_k;
check_simulation_parameters(&sim);
@@ -240,32 +242,32 @@ main(int argc, char* argv[])
t_begin = clock();
#endif
- /* deep deformatin does not occur with approximately zero amplitude in
- * water pressure forcing, or if the stress gradient is positive at
- * zero depth */
- if (fabs(sim.p_f_mod_ampl) > 1e-6 &&
- eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) {
+ /* deep deformatin does not occur with approximately zero
+ * amplitude in water pressure forcing, or if the stress gradient
+ * is positive at zero depth */
+ if (fabs(sim.p_f_mod_ampl) > 1e-6 &&
+ eff_normal_stress_gradient(&sim, d_s, 0.0) < 0.0) {
depth_limit1 = 0.0;
- depth_limit2 = 5.0*d_s;
+ depth_limit2 = 5.0 * d_s;
depth = zbrent(&sim,
- d_s,
- (double (*)(struct simulation*, double, double))
- eff_normal_stress_gradient,
- depth_limit1,
- depth_limit2,
- TOL);
+ d_s,
+ (double (*) (struct simulation *, double, double))
+ eff_normal_stress_gradient,
+ depth_limit1,
+ depth_limit2,
+ TOL);
}
-
#ifdef BENCHMARK_PERFORMANCE
t_end = clock();
- t_elapsed = (double)(t_end - t_begin)/CLOCKS_PER_SEC;
+ t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC;
printf("time spent = %g s\n", t_elapsed);
#endif
printf("%.17g\t%.17g\n", depth, d_s);
free_arrays(&sim);
+
return 0;
}
diff --git a/shear_flux.c b/shear_flux.c
@@ -29,7 +29,7 @@ find_flux(FILE *f)
flux = 0.0;
while (fscanf(f, "%lf\t%lf%*[^\n]", &pos, &vel) == 2) {
if (i++ > 0)
- flux += (vel + vel_prev)/2.0*(pos - pos_prev);
+ flux += (vel + vel_prev) / 2.0 * (pos - pos_prev);
pos_prev = pos;
vel_prev = vel;
}
@@ -40,7 +40,7 @@ find_flux(FILE *f)
}
int
-main(int argc, char* argv[])
+main(int argc, char *argv[])
{
int i;
FILE *fp;
diff --git a/simulation.c b/simulation.c
@@ -35,15 +35,15 @@ init_sim(struct simulation *sim)
sim->v_x_fix = NAN;
sim->v_x_limit = NAN;
- sim->nz = -1; /* cell size equals grain size if negative */
+ sim->nz = -1; /* cell size equals grain size if negative */
- /* lower values of A mean that the velocity curve can have sharper curves,
- * e.g. at the transition from μ ≈ μ_s */
- sim->A = 0.40; /* Loose fit to Damsgaard et al 2013 */
+ /* lower values of A mean that the velocity curve can have sharper
+ * curves, e.g. at the transition from μ ≈ μ_s */
+ sim->A = 0.40; /* Loose fit to Damsgaard et al 2013 */
- /* lower values of b mean larger shear velocity for a given stress ratio
- * above yield */
- sim->b = 0.9377; /* Henann and Kamrin 2016 */
+ /* lower values of b mean larger shear velocity for a given stress
+ * ratio above yield */
+ sim->b = 0.9377; /* Henann and Kamrin 2016 */
/* Henann and Kamrin 2016 */
/* sim->mu_s = 0.3819; */
@@ -53,7 +53,7 @@ init_sim(struct simulation *sim)
sim->mu_s = tan(DEG2RAD(22.0));
sim->C = 0.0;
sim->phi = initval(0.25, 1);
- sim->d = 0.04; /* Damsgaard et al 2013 */
+ sim->d = 0.04; /* Damsgaard et al 2013 */
sim->transient = 0;
sim->phi_min = 0.20;
@@ -79,7 +79,7 @@ init_sim(struct simulation *sim)
/* sim->d = ??; */
/* grain material density [kg/m^3] */
- sim->rho_s = 2.6e3; /* Damsgaard et al 2013 */
+ sim->rho_s = 2.6e3; /* Damsgaard et al 2013 */
/* spatial settings */
sim->origo_z = 0.0;
@@ -102,8 +102,8 @@ init_sim(struct simulation *sim)
/* sim->mu_f = 1.0-3; */
/* Water at 0 deg C */
- sim->beta_f = 3.9e-10; /* doi:10.1063/1.1679903 */
- sim->mu_f = 1.787e-3; /* Cuffey and Paterson 2010 */
+ sim->beta_f = 3.9e-10; /* doi:10.1063/1.1679903 */
+ sim->mu_f = 1.787e-3; /* Cuffey and Paterson 2010 */
/* Damsgaard et al 2015 */
sim->k = initval(1.9e-15, 1);
@@ -131,33 +131,32 @@ prepare_arrays(struct simulation *sim)
{
if (sim->nz < 2) {
fprintf(stderr, "error: grid size (nz) must be at least 2 but is %d\n",
- sim->nz);
+ sim->nz);
exit(1);
}
-
free(sim->phi);
free(sim->k);
- sim->z = linspace(sim->origo_z, /* spatial coordinates */
- sim->origo_z + sim->L_z,
- sim->nz);
- sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */
- sim->mu = zeros(sim->nz); /* stress ratio */
- sim->mu_c = zeros(sim->nz); /* critical-state stress ratio */
- sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */
- sim->sigma_n = zeros(sim->nz); /* normal stess */
- sim->p_f_ghost = zeros(sim->nz+2); /* fluid pressure with ghost nodes */
- sim->phi = zeros(sim->nz); /* porosity */
- sim->phi_c = zeros(sim->nz); /* critical-state porosity */
- sim->phi_dot = zeros(sim->nz); /* rate of porosity change */
- sim->k = zeros(sim->nz); /* permeability */
- sim->xi = zeros(sim->nz); /* cooperativity length */
- sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */
- sim->v_x = zeros(sim->nz); /* shear velocity */
- sim->g_local = zeros(sim->nz); /* local fluidity */
- sim->g_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */
- sim->I = zeros(sim->nz); /* inertia number */
- sim->tan_psi = zeros(sim->nz); /* tan(dilatancy_angle) */
+ sim->z = linspace(sim->origo_z, /* spatial coordinates */
+ sim->origo_z + sim->L_z,
+ sim->nz);
+ sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */
+ sim->mu = zeros(sim->nz); /* stress ratio */
+ sim->mu_c = zeros(sim->nz); /* critical-state stress ratio */
+ sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */
+ sim->sigma_n = zeros(sim->nz); /* normal stess */
+ sim->p_f_ghost = zeros(sim->nz + 2); /* fluid pressure with ghost nodes */
+ sim->phi = zeros(sim->nz); /* porosity */
+ sim->phi_c = zeros(sim->nz); /* critical-state porosity */
+ sim->phi_dot = zeros(sim->nz); /* rate of porosity change */
+ sim->k = zeros(sim->nz);/* permeability */
+ sim->xi = zeros(sim->nz); /* cooperativity length */
+ sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */
+ sim->v_x = zeros(sim->nz); /* shear velocity */
+ sim->g_local = zeros(sim->nz); /* local fluidity */
+ sim->g_ghost = zeros(sim->nz + 2); /* fluidity with ghost nodes */
+ sim->I = zeros(sim->nz);/* inertia number */
+ sim->tan_psi = zeros(sim->nz); /* tan(dilatancy_angle) */
}
void
@@ -184,29 +183,30 @@ free_arrays(struct simulation *sim)
static void
warn_parameter_value(const char message[],
- const double value,
- int* return_status)
+ const double value,
+ int *return_status)
{
fprintf(stderr,
"check_simulation_parameters: %s (%.17g)\n",
- message,
- value);
+ message, value);
*return_status = 1;
}
static void
-check_float(const char name[], const double value, int* return_status)
+check_float(const char name[], const double value, int *return_status)
{
#ifdef SHOW_PARAMETERS
printf("%30s: %.17g\n", name, value);
#endif
if (isnan(value)) {
char message[100];
+
snprintf(message, sizeof(message), "%s is NaN", name);
warn_parameter_value(message, value, return_status);
*return_status = 1;
} else if (isinf(value)) {
char message[100];
+
snprintf(message, sizeof(message), "%s is infinite", name);
warn_parameter_value(message, value, return_status);
*return_status = 1;
@@ -225,14 +225,14 @@ check_simulation_parameters(struct simulation *sim)
check_float("sim->P_wall", sim->P_wall, &return_status);
if (sim->P_wall < 0.0)
warn_parameter_value("sim->P_wall is negative", sim->P_wall,
- &return_status);
+ &return_status);
check_float("sim->v_x_bot", sim->v_x_bot, &return_status);
check_float("sim->mu_wall", sim->mu_wall, &return_status);
if (sim->mu_wall < 0.0)
- warn_parameter_value("sim->mu_wall is negative", sim->mu_wall,
- &return_status);
+ warn_parameter_value("sim->mu_wall is negative", sim->mu_wall,
+ &return_status);
check_float("sim->A", sim->A, &return_status);
if (sim->A < 0.0)
@@ -245,125 +245,124 @@ check_simulation_parameters(struct simulation *sim)
check_float("sim->mu_s", sim->mu_s, &return_status);
if (sim->mu_s < 0.0)
warn_parameter_value("sim->mu_s is negative", sim->mu_s,
- &return_status);
+ &return_status);
check_float("sim->C", sim->C, &return_status);
check_float("sim->d", sim->d, &return_status);
if (sim->d <= 0.0)
warn_parameter_value("sim->d is not a positive number", sim->d,
- &return_status);
+ &return_status);
check_float("sim->rho_s", sim->rho_s, &return_status);
if (sim->rho_s <= 0.0)
warn_parameter_value("sim->rho_s is not a positive number", sim->rho_s,
- &return_status);
+ &return_status);
if (sim->nz <= 0)
warn_parameter_value("sim->nz is not a positive number", sim->nz,
- &return_status);
+ &return_status);
check_float("sim->origo_z", sim->origo_z, &return_status);
check_float("sim->L_z", sim->L_z, &return_status);
if (sim->L_z <= sim->origo_z)
warn_parameter_value("sim->L_z is smaller or equal to sim->origo_z",
- sim->L_z, &return_status);
+ sim->L_z, &return_status);
if (sim->nz <= 0)
warn_parameter_value("sim->nz is not a positive number", sim->nz,
- &return_status);
+ &return_status);
check_float("sim->dz", sim->dz, &return_status);
if (sim->dz <= 0.0)
warn_parameter_value("sim->dz is not a positive number", sim->dz,
- &return_status);
+ &return_status);
check_float("sim->t", sim->t, &return_status);
if (sim->t < 0.0)
warn_parameter_value("sim->t is a negative number",
- sim->t, &return_status);
+ sim->t, &return_status);
check_float("sim->t_end", sim->t_end, &return_status);
if (sim->t > sim->t_end)
warn_parameter_value("sim->t_end is smaller than sim->t",
- sim->t, &return_status);
+ sim->t, &return_status);
check_float("sim->dt", sim->dt, &return_status);
if (sim->dt <= 0.0)
warn_parameter_value("sim->dt is not a positive number",
- sim->dt, &return_status);
+ sim->dt, &return_status);
check_float("sim->file_dt", sim->file_dt, &return_status);
if (sim->file_dt < 0.0)
warn_parameter_value("sim->file_dt is a negative number",
- sim->file_dt, &return_status);
+ sim->file_dt, &return_status);
check_float("sim->phi[0]", sim->phi[0], &return_status);
if (sim->phi[0] < 0.0 || sim->phi[0] > 1.0)
warn_parameter_value("sim->phi[0] is not within [0;1]",
- sim->phi[0], &return_status);
+ sim->phi[0], &return_status);
check_float("sim->phi_min", sim->phi_min, &return_status);
if (sim->phi_min < 0.0 || sim->phi_min > 1.0)
warn_parameter_value("sim->phi_min is not within [0;1]",
- sim->phi_min, &return_status);
+ sim->phi_min, &return_status);
check_float("sim->phi_max", sim->phi_max, &return_status);
if (sim->phi_max < 0.0 || sim->phi_max > 1.0)
warn_parameter_value("sim->phi_max is not within [0;1]",
- sim->phi_max, &return_status);
+ sim->phi_max, &return_status);
check_float("sim->dilatancy_angle", sim->dilatancy_angle, &return_status);
if (sim->dilatancy_angle < 0.0 || sim->dilatancy_angle > 100.0)
warn_parameter_value("sim->dilatancy_angle is not within [0;100]",
- sim->dilatancy_angle, &return_status);
+ sim->dilatancy_angle, &return_status);
if (sim->fluid != 0 && sim->fluid != 1)
warn_parameter_value("sim->fluid has an invalid value",
- (double)sim->fluid, &return_status);
+ (double) sim->fluid, &return_status);
if (sim->transient != 0 && sim->transient != 1)
warn_parameter_value("sim->transient has an invalid value",
- (double)sim->transient, &return_status);
+ (double) sim->transient, &return_status);
if (sim->fluid) {
check_float("sim->p_f_mod_ampl", sim->p_f_mod_ampl, &return_status);
if (sim->p_f_mod_ampl < 0.0)
warn_parameter_value("sim->p_f_mod_ampl is not a zero or positive",
- sim->p_f_mod_ampl, &return_status);
+ sim->p_f_mod_ampl, &return_status);
if (sim->P_wall - sim->p_f_mod_ampl < 0.0)
warn_parameter_value("sim->P_wall - sim->p_f_mod_ampl is negative",
- sim->P_wall - sim->p_f_mod_ampl,
- &return_status);
+ sim->P_wall - sim->p_f_mod_ampl,
+ &return_status);
check_float("sim->p_f_mod_freq", sim->p_f_mod_freq, &return_status);
- if (sim->p_f_mod_freq < 0.0)
- warn_parameter_value("sim->p_f_mod_freq is not a zero or positive",
- sim->p_f_mod_freq, &return_status);
+ if (sim->p_f_mod_freq < 0.0)
+ warn_parameter_value("sim->p_f_mod_freq is not a zero or positive",
+ sim->p_f_mod_freq, &return_status);
check_float("sim->beta_f", sim->beta_f, &return_status);
- if (sim->beta_f <= 0.0)
- warn_parameter_value("sim->beta_f is not positive",
- sim->beta_f, &return_status);
+ if (sim->beta_f <= 0.0)
+ warn_parameter_value("sim->beta_f is not positive",
+ sim->beta_f, &return_status);
check_float("sim->mu_f", sim->mu_f, &return_status);
- if (sim->mu_f <= 0.0)
- warn_parameter_value("sim->mu_f is not positive",
- sim->mu_f, &return_status);
+ if (sim->mu_f <= 0.0)
+ warn_parameter_value("sim->mu_f is not positive",
+ sim->mu_f, &return_status);
check_float("sim->rho_f", sim->rho_f, &return_status);
- if (sim->rho_f <= 0.0)
- warn_parameter_value("sim->rho_f is not positive",
- sim->rho_f, &return_status);
+ if (sim->rho_f <= 0.0)
+ warn_parameter_value("sim->rho_f is not positive",
+ sim->rho_f, &return_status);
check_float("sim->k[0]", sim->k[0], &return_status);
- if (sim->k[0] <= 0.0)
- warn_parameter_value("sim->k[0] is not positive",
- sim->k[0], &return_status);
+ if (sim->k[0] <= 0.0)
+ warn_parameter_value("sim->k[0] is not positive",
+ sim->k[0], &return_status);
}
-
if (return_status != 0) {
fprintf(stderr, "error: aborting due to invalid parameter choices\n");
exit(return_status);
@@ -374,10 +373,11 @@ void
lithostatic_pressure_distribution(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
+
+ for (i = 0; i < sim->nz; ++i)
sim->sigma_n[i] = sim->P_wall +
- (1.0 - sim->phi[i])*sim->rho_s*sim->G*
- (sim->L_z - sim->z[i]);
+ (1.0 - sim->phi[i]) * sim->rho_s * sim->G *
+ (sim->L_z - sim->z[i]);
}
/* NEW FUNCTIONS START */
@@ -386,42 +386,47 @@ void
compute_inertia_number(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->I[i] = sim->gamma_dot_p[i]*sim->d
- /sqrt(sim->sigma_n_eff[i]/sim->rho_s);
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->I[i] = sim->gamma_dot_p[i] * sim->d
+ / sqrt(sim->sigma_n_eff[i] / sim->rho_s);
}
+
void
compute_critical_state_porosity(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min)*sim->I[i];
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min) *sim->I[i];
}
void
compute_critical_state_friction(struct simulation *sim)
{
int i;
+
if (sim->fluid)
- for (i=0; i<sim->nz; ++i)
- sim->mu_c[i] = sim->mu_wall/
- (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top));
+ for (i = 0; i < sim->nz; ++i)
+ sim->mu_c[i] = sim->mu_wall /
+ (sim->sigma_n_eff[i] / (sim->P_wall - sim->p_f_top));
else
- for (i=0; i<sim->nz; ++i)
- sim->mu_c[i] = sim->mu_wall/
- (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G*
- (sim->L_z - sim->z[i])/sim->P_wall);
+ for (i = 0; i < sim->nz; ++i)
+ sim->mu_c[i] = sim->mu_wall /
+ (1.0 + (1.0 - sim->phi[i]) * sim->rho_s * sim->G *
+ (sim->L_z - sim->z[i]) / sim->P_wall);
}
void
compute_friction(struct simulation *sim)
{
int i;
+
if (sim->transient)
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
sim->mu[i] = sim->tan_psi[i] + sim->mu_c[i];
else
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
sim->mu[i] = sim->mu_c[i];
}
@@ -429,17 +434,19 @@ void
compute_tan_dilatancy_angle(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->tan_psi[i] = sim->dilatancy_angle*(sim->phi_c[i] - sim->phi[i]);
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->tan_psi[i] = sim->dilatancy_angle * (sim->phi_c[i] - sim->phi[i]);
}
void
compute_porosity_change(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i) {
- sim->phi_dot[i] = sim->tan_psi[i]*sim->gamma_dot_p[i]*sim->phi[i];
- sim->phi[i] += sim->phi_dot[i]*sim->dt;
+
+ for (i = 0; i < sim->nz; ++i) {
+ sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i];
+ sim->phi[i] += sim->phi_dot[i] * sim->dt;
}
}
@@ -447,42 +454,29 @@ void
compute_permeability(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->k[i] = pow(sim->d, 2.0)/180.0
- *pow(sim->phi[i], 3.0)
- /pow(1.0 - sim->phi[i], 2.0);
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->k[i] = pow(sim->d, 2.0) / 180.0
+ * pow(sim->phi[i], 3.0)
+ / pow(1.0 - sim->phi[i], 2.0);
}
/* NEW FUNCTIONS END */
-/*void
-compute_friction(struct simulation *sim)
-{
- int i;
- if (sim->fluid)
- for (i=0; i<sim->nz; ++i)
- sim->mu[i] = sim->mu_wall/
- (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top));
- else
- for (i=0; i<sim->nz; ++i)
- sim->mu[i] = sim->mu_wall/
- (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G*
- (sim->L_z - sim->z[i])/sim->P_wall);
-}*/
-
static double
shear_strain_rate_plastic(const double fluidity, const double friction)
{
- return fluidity*friction;
+ return fluidity * friction;
}
void
compute_shear_strain_rate_plastic(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
- sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[i+1],
- sim->mu[i]);
+
+ for (i = 0; i < sim->nz; ++i)
+ sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[i + 1],
+ sim->mu[i]);
}
void
@@ -494,19 +488,20 @@ compute_shear_velocity(struct simulation *sim)
/* Dirichlet BC at bottom */
sim->v_x[0] = sim->v_x_bot;
- for (i=1; i<sim->nz; ++i)
- sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz;
+ for (i = 1; i < sim->nz; ++i)
+ sim->v_x[i] = sim->v_x[i - 1] + sim->gamma_dot_p[i] * sim->dz;
}
void
compute_effective_stress(struct simulation *sim)
{
int i;
+
if (sim->fluid)
- for (i=0; i<sim->nz; ++i)
- sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i+1];
+ for (i = 0; i < sim->nz; ++i)
+ sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i + 1];
else
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
sim->sigma_n_eff[i] = sim->sigma_n[i];
}
@@ -518,14 +513,15 @@ cooperativity_length(const double A,
const double mu_s,
const double C)
{
- return A*d/sqrt(fabs((mu - C/p) - mu_s));
+ return A * d / sqrt(fabs((mu - C / p) - mu_s));
}
void
compute_cooperativity_length(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
+
+ for (i = 0; i < sim->nz; ++i)
sim->xi[i] = cooperativity_length(sim->A,
sim->d,
sim->mu[i],
@@ -543,17 +539,18 @@ local_fluidity(const double p,
const double rho_s,
const double d)
{
- if (mu - C/p <= mu_s)
- return 0.0;
+ if (mu - C / p <= mu_s)
+ return 0.0;
else
- return sqrt(p/rho_s*d*d)*((mu - C/p) - mu_s)/(b*mu);
+ return sqrt(p / rho_s * d * d) * ((mu - C / p) - mu_s) / (b * mu);
}
void
compute_local_fluidity(struct simulation *sim)
{
int i;
- for (i=0; i<sim->nz; ++i)
+
+ for (i = 0; i < sim->nz; ++i)
sim->g_local[i] = local_fluidity(sim->sigma_n_eff[i],
sim->mu[i],
sim->mu_s,
@@ -564,16 +561,16 @@ compute_local_fluidity(struct simulation *sim)
}
void
-set_bc_neumann(double* a,
+set_bc_neumann(double *a,
const int nz,
const int boundary,
const double df,
const double dx)
{
if (boundary == -1)
- a[0] = a[1] + df*dx;
+ a[0] = a[1] + df * dx;
else if (boundary == +1)
- a[nz+1] = a[nz] - df*dx;
+ a[nz + 1] = a[nz] - df * dx;
else {
fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary);
exit(1);
@@ -581,7 +578,7 @@ set_bc_neumann(double* a,
}
void
-set_bc_dirichlet(double* a,
+set_bc_dirichlet(double *a,
const int nz,
const int boundary,
const double value)
@@ -589,7 +586,7 @@ set_bc_dirichlet(double* a,
if (boundary == -1)
a[0] = value;
else if (boundary == +1)
- a[nz+1] = value;
+ a[nz + 1] = value;
else {
fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary);
exit(1);
@@ -599,33 +596,33 @@ set_bc_dirichlet(double* a,
double
residual_normalized(double new, double old)
{
- return fabs((new - old)/(old + 1e-16));
+ return fabs((new - old) / (old + 1e-16));
}
static void
poisson_solver_1d_cell_update(int i,
- const double* g_in,
- const double* g_local,
- double* g_out,
- double* r_norm,
+ const double *g_in,
+ const double *g_local,
+ double *g_out,
+ double *r_norm,
const double dz,
- const double* xi)
+ const double *xi)
{
- double coorp_term = dz*dz/(2.0*pow(xi[i], 2.0));
+ double coorp_term = dz * dz / (2.0 * pow(xi[i], 2.0));
- g_out[i+1] = 1.0/(1.0 + coorp_term)
- *(coorp_term*g_local[i] + g_in[i+2]/2.0 + g_in[i]/2.0);
+ g_out[i + 1] = 1.0 / (1.0 + coorp_term)
+ * (coorp_term * g_local[i] + g_in[i + 2] / 2.0 + g_in[i] / 2.0);
/* TODO: try out residual_normalized(..) */
- r_norm[i] = pow(g_out[i+1] - g_in[i+1], 2.0)
- /(pow(g_out[i+1], 2.0) + 1e-16);
+ r_norm[i] = pow(g_out[i + 1] - g_in[i + 1], 2.0)
+ / (pow(g_out[i + 1], 2.0) + 1e-16);
#ifdef DEBUG
printf("-- %d --------------\n", i);
printf("coorp_term: %g\n", coorp_term);
printf(" g_local: %g\n", g_local[i]);
- printf(" g_in: %g\n", g_in[i+1]);
- printf(" g_out: %g\n", g_out[i+1]);
+ printf(" g_in: %g\n", g_in[i + 1]);
+ printf(" g_out: %g\n", g_out[i + 1]);
printf(" r_norm: %g\n", r_norm[i]);
#endif
}
@@ -637,9 +634,9 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
{
int iter, i;
double r_norm_max = NAN;
- double *g_ghost_out = empty(sim->nz+2), *r_norm = empty(sim->nz);
+ double *g_ghost_out = empty(sim->nz + 2), *r_norm = empty(sim->nz);
- for (iter=0; iter<max_iter; ++iter) {
+ for (iter = 0; iter < max_iter; ++iter) {
#ifdef DEBUG
printf("\n@@@ ITERATION %d @@@\n", iter);
#endif
@@ -650,7 +647,7 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
/* Neumann BCs resemble free surfaces */
/* set_bc_neumann(sim->g_ghost, sim->nz, +1, 0.0, sim->dz); */
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
poisson_solver_1d_cell_update(i,
sim->g_ghost,
sim->g_local,
@@ -660,14 +657,15 @@ implicit_1d_jacobian_poisson_solver(struct simulation *sim,
sim->xi);
r_norm_max = max(r_norm, sim->nz);
- copy_values(g_ghost_out, sim->g_ghost, sim->nz+2);
+ copy_values(g_ghost_out, sim->g_ghost, sim->nz + 2);
if (r_norm_max <= rel_tol) {
set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0);
free(g_ghost_out);
free(r_norm);
- /* printf(".. Solution converged after %d iterations\n", iter); */
+ /* printf(".. Solution converged after %d
+ * iterations\n", iter); */
return 0;
}
}
@@ -687,7 +685,7 @@ write_output_file(struct simulation *sim, const int normalize)
FILE *fp;
snprintf(outfile, sizeof(outfile), "%s.output%05d.txt",
- sim->name, sim->n_file++);
+ sim->name, sim->n_file++);
if ((fp = fopen(outfile, "w")) != NULL) {
print_output(sim, fp, normalize);
@@ -699,7 +697,7 @@ write_output_file(struct simulation *sim, const int normalize)
}
void
-print_output(struct simulation *sim, FILE* fp, const int norm)
+print_output(struct simulation *sim, FILE *fp, const int norm)
{
int i;
double *v_x_out;
@@ -709,7 +707,7 @@ print_output(struct simulation *sim, FILE* fp, const int norm)
else
v_x_out = copy(sim->v_x, sim->nz);
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
fprintf(fp,
"%.17g\t%.17g\t%.17g\t"
"%.17g\t%.17g\t%.17g\t"
@@ -717,20 +715,20 @@ print_output(struct simulation *sim, FILE* fp, const int norm)
sim->z[i],
v_x_out[i],
sim->sigma_n_eff[i],
- sim->p_f_ghost[i+1],
+ sim->p_f_ghost[i + 1],
sim->mu[i],
sim->gamma_dot_p[i],
sim->phi[i],
sim->I[i],
- sim->mu[i]*sim->sigma_n_eff[i]);
+ sim->mu[i] * sim->sigma_n_eff[i]);
free(v_x_out);
}
int
coupled_shear_solver(struct simulation *sim,
- const int max_iter,
- const double rel_tol)
+ const int max_iter,
+ const double rel_tol)
{
int i, coupled_iter, stress_iter = 0;
double *r_norm, r_norm_max, *oldval;
@@ -740,30 +738,29 @@ coupled_shear_solver(struct simulation *sim,
r_norm = empty(sim->nz);
oldval = empty(sim->nz);
}
-
- do { /* stress iterations */
+ do { /* stress iterations */
coupled_iter = 0;
- do { /* coupled iterations */
+ do { /* coupled iterations */
if (sim->transient) {
copy_values(sim->gamma_dot_p, oldval, sim->nz);
/* step 1 */
- compute_inertia_number(sim); /* Eq. 1 */
+ compute_inertia_number(sim); /* Eq. 1 */
/* step 2 */
compute_critical_state_porosity(sim); /* Eq. 2 */
/* step 3 */
- compute_tan_dilatancy_angle(sim); /* Eq. 5 */
+ compute_tan_dilatancy_angle(sim); /* Eq. 5 */
}
- compute_critical_state_friction(sim); /* Eq. 7 */
+ compute_critical_state_friction(sim); /* Eq. 7 */
/* step 4 */
if (sim->transient) {
- compute_porosity_change(sim); /* Eq. 3 */
- compute_permeability(sim); /* Eq. 6 */
+ compute_porosity_change(sim); /* Eq. 3 */
+ compute_permeability(sim); /* Eq. 6 */
}
- compute_friction(sim); /* Eq. 4 */
+ compute_friction(sim); /* Eq. 4 */
/* step 5, Eq. 13 */
if (sim->fluid)
@@ -771,16 +768,16 @@ coupled_shear_solver(struct simulation *sim,
exit(11);
/* step 6 */
- compute_effective_stress(sim); /* Eq. 9 */
+ compute_effective_stress(sim); /* Eq. 9 */
/* step 7 */
- compute_local_fluidity(sim); /* Eq. 10 */
- compute_cooperativity_length(sim); /* Eq. 12 */
+ compute_local_fluidity(sim); /* Eq. 10 */
+ compute_cooperativity_length(sim); /* Eq. 12 */
/* step 8, Eq. 11 */
if (implicit_1d_jacobian_poisson_solver(sim,
- MAX_ITER_GRANULAR,
- RTOL_GRANULAR))
+ MAX_ITER_GRANULAR,
+ RTOL_GRANULAR))
exit(12);
/* step 9 */
@@ -789,7 +786,7 @@ coupled_shear_solver(struct simulation *sim,
compute_shear_velocity(sim);
#ifdef DEBUG
- for (i=0; i<sim->nz) {
+ for (i = 0; i < sim->nz) {
printf("\nsim->t = %g s\n", sim->t);
printf("sim->I[%d] = %g\n", i, sim->I[i]);
printf("sim->phi_c[%d] = %g\n", i, sim->phi_c[i]);
@@ -800,9 +797,10 @@ coupled_shear_solver(struct simulation *sim,
printf("sim->mu[%d] = %g\n", i, sim->mu[i]);
}
#endif
- /* stable velocity field == coupled solution converged */
+ /* stable velocity field == coupled solution
+ * converged */
if (sim->transient) {
- for (i=0; i<sim->nz; ++i)
+ for (i = 0; i < sim->nz; ++i)
r_norm[i] = residual_normalized(sim->gamma_dot_p[i], oldval[i]);
r_norm_max = max(r_norm, sim->nz);
if (r_norm_max <= rel_tol)
@@ -813,9 +811,9 @@ coupled_shear_solver(struct simulation *sim,
free(oldval);
fprintf(stderr, "coupled_shear_solver: ");
fprintf(stderr, "Transient solution did not converge "
- "after %d iterations\n", coupled_iter);
+ "after %d iterations\n", coupled_iter);
fprintf(stderr, ".. Residual normalized error: %g\n",
- r_norm_max);
+ r_norm_max);
return 1;
}
}
@@ -823,28 +821,26 @@ coupled_shear_solver(struct simulation *sim,
if (!isnan(sim->v_x_limit) || !isnan(sim->v_x_fix)) {
if (!isnan(sim->v_x_limit)) {
- vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz-1])
- /(sim->v_x[sim->nz-1] + 1e-12);
+ vel_res_norm = (sim->v_x_limit - sim->v_x[sim->nz - 1])
+ / (sim->v_x[sim->nz - 1] + 1e-12);
if (vel_res_norm > 0.0)
vel_res_norm = 0.0;
} else {
- vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz-1])
- /(sim->v_x[sim->nz-1] + 1e-12);
+ vel_res_norm = (sim->v_x_fix - sim->v_x[sim->nz - 1])
+ / (sim->v_x[sim->nz - 1] + 1e-12);
}
- sim->mu_wall *= 1.0 + (vel_res_norm*1e-2);
+ sim->mu_wall *= 1.0 + (vel_res_norm * 1e-2);
}
-
if (++stress_iter > MAX_ITER_STRESS) {
fprintf(stderr, "error: stress solution did not converge:\n");
fprintf(stderr, "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
- "vel_res_norm=%g, mu_wall=%g\n",
- sim->v_x[sim->nz-1], sim->v_x_fix, sim->v_x_limit,
- vel_res_norm, sim->mu_wall);
+ "vel_res_norm=%g, mu_wall=%g\n",
+ sim->v_x[sim->nz - 1], sim->v_x_fix, sim->v_x_limit,
+ vel_res_norm, sim->mu_wall);
return 10;
}
-
} while ((!isnan(sim->v_x_fix) || !isnan(sim->v_x_limit))
- && fabs(vel_res_norm) > RTOL_STRESS);
+ && fabs(vel_res_norm) > RTOL_STRESS);
if (!isnan(sim->v_x_limit))
sim->mu_wall = mu_wall_orig;
@@ -853,18 +849,17 @@ coupled_shear_solver(struct simulation *sim,
free(r_norm);
free(oldval);
}
-
return 0;
}
double
-find_flux(const struct simulation* sim)
+find_flux(const struct simulation *sim)
{
int i;
double flux = 0.0;
- for (i=1; i<sim->nz; ++i)
- flux += (sim->v_x[i] + sim->v_x[i-1])/2.0*sim->dz;
+ for (i = 1; i < sim->nz; ++i)
+ flux += (sim->v_x[i] + sim->v_x[i - 1]) / 2.0 * sim->dz;
return flux;
}
