tUse tabs for indentation and spaces for alignment - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear) (HTM) git clone git://src.adamsgaard.dk/cngf-pf (DIR) Log (DIR) Files (DIR) Refs (DIR) README (DIR) LICENSE --- (DIR) commit 4a6213835f491e0753ed1a3f864380a75a0d30cd (DIR) parent 8e120cffa8aad6b9eab4ff77a4dfdc0ab5402dd2 (HTM) Author: Anders Damsgaard <anders@adamsgaard.dk> Date: Tue, 25 Jun 2019 16:41:55 +0200 Use tabs for indentation and spaces for alignment Diffstat: M arrays.c | 148 ++++++++++++++++---------------- M fluid.c | 281 +++++++++++++++---------------- M main.c | 553 ++++++++++++++++--------------- M simulation.c | 683 +++++++++++++++---------------- 4 files changed, 827 insertions(+), 838 deletions(-) --- (DIR) diff --git a/arrays.c b/arrays.c t@@ -6,158 +6,158 @@ /* 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, const unsigned int k, - const unsigned int nx, const unsigned int ny) +unsigned int idx3(const unsigned int i, + const unsigned int j, + const unsigned int k, + 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 */ -unsigned int idx3g( - const unsigned int i, const unsigned int j, const unsigned int k, - const unsigned int nx, const unsigned int ny) +unsigned int idx3g(const unsigned int i, + const unsigned int j, + const unsigned int k, + 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 */ -unsigned int idx2( - const unsigned int i, const unsigned int j, const unsigned int nx) +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 */ -unsigned int idx2g( - const unsigned int i, const unsigned int j, const unsigned int nx) +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 */ 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* linspace(const double lower, const double upper, const int n) { - double *x = malloc(n*sizeof(double)); - double dx = (upper - lower)/(double)(n-1); - for (int i=0; i<n; ++i) - x[i] = lower + dx*i; - return x; + double *x = malloc(n*sizeof(double)); + double dx = (upper - lower)/(double)(n-1); + for (int i=0; i<n; ++i) + x[i] = lower + dx*i; + return x; } /* Return an array of `n` values with the value 0.0 */ double* zeros(const int n) { - double *x = malloc(n*sizeof(double)); - for (int i=0; i<n; ++i) - x[i] = 0.0; - return x; + double *x = malloc(n*sizeof(double)); + for (int i=0; i<n; ++i) + x[i] = 0.0; + return x; } /* Return an array of `n` values with the value 1.0 */ double* ones(const int n) { - double *x = malloc(n*sizeof(double)); - for (int i=0; i<n; ++i) - x[i] = 1.0; - return x; + double *x = malloc(n*sizeof(double)); + for (int i=0; i<n; ++i) + x[i] = 1.0; + return x; } /* Return an array of `n` values with a specified value */ double* initval(const double value, const int n) { - double *x = malloc(n*sizeof(double)); - for (int i=0; i<n; ++i) - x[i] = value; - return x; + double *x = malloc(n*sizeof(double)); + for (int i=0; i<n; ++i) + x[i] = value; + return x; } /* Return an array of `n` uninitialized values */ double* empty(const int 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) { - double maxval = -INFINITY; - for (int i=0; i<n; ++i) - if (a[i] > maxval) - maxval = a[i]; - return maxval; + double maxval = -INFINITY; + for (int 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) { - double minval = +INFINITY; - for (int i=0; i<n; ++i) - if (a[i] < minval) - minval = a[i]; - return minval; + double minval = +INFINITY; + for (int i=0; i<n; ++i) + if (a[i] < minval) + minval = a[i]; + return minval; } void print_array(const double* a, const int n) { - for (int i=0; i<n; ++i) - printf("%.17g\n", a[i]); + for (int i=0; i<n; ++i) + printf("%.17g\n", a[i]); } void print_arrays(const double* a, const double* b, const int n) { - for (int i=0; i<n; ++i) - printf("%.17g\t%.17g\n", a[i], b[i]); + for (int 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) { - double max_b = max(b, n); - for (int i=0; i<n; ++i) - printf("%.17g\t%.17g\n", a[i], b[i]/max_b); + double max_b = max(b, n); + for (int 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, - const int n) +void print_three_arrays(const double* a, + const double* b, + const double* c, + const int n) { - for (int i=0; i<n; ++i) - printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); + for (int 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) +void fprint_arrays(FILE* fp, const double* a, const double* b, const int n) { - for (int i=0; i<n; ++i) - fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]); + for (int 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, - const int n) +void fprint_three_arrays(FILE* fp, + const double* a, + const double* b, + const double* c, + const int n) { - for (int i=0; i<n; ++i) - fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); + for (int 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) { - for (int i=0; i<n; ++i) - out[i] = in[i]; + for (int i=0; i<n; ++i) + out[i] = in[i]; } (DIR) diff --git a/fluid.c b/fluid.c t@@ -5,180 +5,177 @@ void hydrostatic_fluid_pressure_distribution(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) - sim->p_f_ghost[idx1g(i)] = sim->p_f_top + - sim->phi[i]*sim->rho_f*sim->G*(sim->L_z - sim->z[i]); + for (int i=0; i<sim->nz; ++i) + sim->p_f_ghost[idx1g(i)] = sim->p_f_top + + sim->phi[i]*sim->rho_f*sim->G* + (sim->L_z - sim->z[i]); } -static double sine_wave( - const double time, - const double amplitude, - const double frequency, - const double phase, - const double base_value) +static double sine_wave(const double time, + const double amplitude, + const double frequency, + const double phase, + const double base_value) { - return amplitude*sin(2.0*PI*frequency*time + phase) + base_value; + return amplitude*sin(2.0*PI*frequency*time + phase) + base_value; } -static double darcy_pressure_change_1d( - const int i, - const int nz, - const double* p_f_ghost_in, - const double* phi, - const double* k, - const double dz, - const double dt, - const double beta_f, - const double mu_f) +static double darcy_pressure_change_1d(const int i, + const int nz, + const double* p_f_ghost_in, + const double* phi, + const double* k, + const double dz, + const double dt, + const double beta_f, + const double mu_f) { - const double p = p_f_ghost_in[idx1g(i)]; - const double p_zn = p_f_ghost_in[idx1g(i-1)]; - const double p_zp = p_f_ghost_in[idx1g(i+1)]; - - const double k_ = k[i]; - double 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]; + const double p = p_f_ghost_in[idx1g(i)]; + const double p_zn = p_f_ghost_in[idx1g(i-1)]; + const double p_zp = p_f_ghost_in[idx1g(i+1)]; + + const double k_ = k[i]; + double 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]; #ifdef DEBUG - printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n", - i, idx1g(i), - p_zn, p, p_zp, - k_zn, k_, k_zp); + printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n", + i, idx1g(i), + p_zn, p, p_zp, + k_zn, k_, k_zp); #endif - const double div_k_grad_p = - (2.0*k_zp*k_/(k_zp + k_) * (p_zp - p)/dz - - 2.0*k_zn*k_/(k_zn + k_) * (p - p_zn)/dz - )/dz; + const double div_k_grad_p = + (2.0*k_zp*k_/(k_zp + k_) * (p_zp - p)/dz - + 2.0*k_zn*k_/(k_zn + k_) * (p - p_zn)/dz + )/dz; #ifdef DEBUG - printf("phi[%d]=%g\tdiv_k_grad_p[%d]=%g\n", - i, phi[i], i, div_k_grad_p); + printf("phi[%d]=%g\tdiv_k_grad_p[%d]=%g\n", + i, phi[i], i, div_k_grad_p); #endif - /* return delta p */ - return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p; + /* return delta p */ + return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p; } -int darcy_solver_1d( - struct simulation* sim, - const int max_iter, - const double rel_tol) +int darcy_solver_1d(struct simulation* sim, + const int max_iter, + const double rel_tol) { - /* compute explicit solution to pressure change */ - double* dp_f_expl = zeros(sim->nz); - for (int i=0; i<sim->nz; ++i) - dp_f_expl[i] = darcy_pressure_change_1d( - i, - sim->nz, - sim->p_f_ghost, - sim->phi, - sim->k, - sim->dz, - sim->dt, - sim->beta_f, - sim->mu_f); - - /* choose integration method, parameter in [0.0; 1.0] - * epsilon = 0.0: explicit - * epsilon = 0.5: Crank-Nicolson - * epsilon = 1.0: implicit */ - const double 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 */ - const double theta = 0.05; - /* const double theta = 1.7; */ - - double p_f; - - /* compute implicit solution to pressure change */ - int iter; - double* dp_f_impl = zeros(sim->nz); - double* p_f_ghost_out = zeros(sim->nz+2); - double* r_norm = zeros(sim->nz); - double r_norm_max = NAN; - double p_f_top = sine_wave( - sim->t, - sim->p_f_mod_ampl, - sim->p_f_mod_freq, - sim->p_f_mod_phase, - sim->p_f_top); - - for (iter=0; iter<max_iter; ++iter) { - - set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); - sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* Include top node in BC */ - set_bc_neumann(sim->p_f_ghost, sim->nz, -1); + /* compute explicit solution to pressure change */ + double* dp_f_expl = zeros(sim->nz); + for (int i=0; i<sim->nz; ++i) + dp_f_expl[i] = darcy_pressure_change_1d(i, + sim->nz, + sim->p_f_ghost, + sim->phi, + sim->k, + sim->dz, + sim->dt, + sim->beta_f, + sim->mu_f); + + /* choose integration method, parameter in [0.0; 1.0] + * epsilon = 0.0: explicit + * epsilon = 0.5: Crank-Nicolson + * epsilon = 1.0: implicit */ + const double 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 */ + const double theta = 0.05; + /* const double theta = 1.7; */ + + double p_f; + + /* compute implicit solution to pressure change */ + int iter; + double* dp_f_impl = zeros(sim->nz); + double* p_f_ghost_out = zeros(sim->nz+2); + double* r_norm = zeros(sim->nz); + double r_norm_max = NAN; + double p_f_top = sine_wave( + sim->t, + sim->p_f_mod_ampl, + sim->p_f_mod_freq, + sim->p_f_mod_phase, + sim->p_f_top); + + for (iter=0; iter<max_iter; ++iter) { + + set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); + sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* Include top node in BC */ + set_bc_neumann(sim->p_f_ghost, sim->nz, -1); #ifdef DEBUG - puts(".. p_f_ghost after BC:"); print_array(sim->p_f_ghost, sim->nz+2); + puts(".. p_f_ghost after BC:"); print_array(sim->p_f_ghost, sim->nz+2); #endif - /* for (int i=0; i<sim->nz; ++i) */ - for (int i=0; i<sim->nz-1; ++i) - dp_f_impl[i] = darcy_pressure_change_1d( - i, - sim->nz, - sim->p_f_ghost, - sim->phi, - sim->k, - sim->dz, - sim->dt, - sim->beta_f, - sim->mu_f); - /* for (int i=0; i<sim->nz; ++i) { */ - for (int i=0; i<sim->nz-1; ++i) { + /* for (int i=0; i<sim->nz; ++i) */ + for (int i=0; i<sim->nz-1; ++i) + dp_f_impl[i] = darcy_pressure_change_1d(i, + sim->nz, + sim->p_f_ghost, + sim->phi, + sim->k, + sim->dz, + sim->dt, + sim->beta_f, + sim->mu_f); + /* for (int i=0; i<sim->nz; ++i) { */ + for (int i=0; i<sim->nz-1; ++i) { #ifdef DEBUG - printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n", - i, dp_f_expl[i], i, dp_f_impl[i]); + printf("dp_f_expl[%d] = %g\ndp_f_impl[%d] = %g\n", + i, dp_f_expl[i], i, dp_f_impl[i]); #endif - p_f = sim->p_f_ghost[idx1g(i)]; + p_f = sim->p_f_ghost[idx1g(i)]; - p_f_ghost_out[idx1g(i)] = p_f - + (1.0 - epsilon)*dp_f_expl[i] + epsilon*dp_f_impl[i]; + p_f_ghost_out[idx1g(i)] = p_f + + (1.0 - epsilon)*dp_f_expl[i] + + epsilon*dp_f_impl[i]; - /* apply relaxation */ - p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta) - + p_f_ghost_out[idx1g(i)]*theta; + /* apply relaxation */ + p_f_ghost_out[idx1g(i)] = p_f*(1.0 - theta) + + p_f_ghost_out[idx1g(i)]*theta; - r_norm[i] = (p_f_ghost_out[idx1g(i)] - p_f)/(p_f + 1e-16); - } + r_norm[i] = (p_f_ghost_out[idx1g(i)] - p_f)/(p_f + 1e-16); + } - r_norm_max = max(r_norm, sim->nz); + r_norm_max = max(r_norm, sim->nz); #ifdef DEBUG - puts(".. p_f_ghost_out:"); print_array(p_f_ghost_out, sim->nz+2); + puts(".. p_f_ghost_out:"); print_array(p_f_ghost_out, sim->nz+2); #endif - copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2); + copy_values(p_f_ghost_out, sim->p_f_ghost, sim->nz+2); #ifdef DEBUG - puts(".. p_f_ghost after update:"); - print_array(sim->p_f_ghost, sim->nz+2); + puts(".. p_f_ghost after update:"); + print_array(sim->p_f_ghost, sim->nz+2); #endif - if (r_norm_max <= rel_tol) { - set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); - sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* top node in BC */ - set_bc_neumann(sim->p_f_ghost, sim->nz, -1); - free(dp_f_expl); - free(dp_f_impl); - free(p_f_ghost_out); - free(r_norm); + if (r_norm_max <= rel_tol) { + set_bc_dirichlet(sim->p_f_ghost, sim->nz, +1, p_f_top); + sim->p_f_ghost[idx1g(sim->nz-1)] = p_f_top; /* top node in BC */ + set_bc_neumann(sim->p_f_ghost, sim->nz, -1); + free(dp_f_expl); + free(dp_f_impl); + free(p_f_ghost_out); + free(r_norm); #ifdef DEBUG - printf(".. Solution converged after %d iterations\n", iter); + printf(".. Solution converged after %d iterations\n", iter); #endif - return 0; - } - } - - free(dp_f_expl); - free(dp_f_impl); - free(p_f_ghost_out); - free(r_norm); - fprintf(stderr, "darcy_solver_1d: "); - fprintf(stderr, "Solution did not converge after %d iterations\n", iter); - fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); - return 1; + return 0; + } + } + + free(dp_f_expl); + free(dp_f_impl); + free(p_f_ghost_out); + free(r_norm); + fprintf(stderr, "darcy_solver_1d: "); + fprintf(stderr, "Solution did not converge after %d iterations\n", iter); + fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); + return 1; } (DIR) diff --git a/main.c b/main.c t@@ -13,317 +13,318 @@ static void usage(void) { - struct simulation sim = init_sim(); - printf("%s: %s [OPTIONS] [NAME]\n" - "runs a simulation and outputs state in files prefixed with NAME.\n" - "If NAME is not specified, the default value '%s' is used.\n" - "optional arguments:\n" - " -N, --normalize normalize output velocity\n" - " -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n" - " -P, --normal-stress VAL normal stress on top [Pa] (default %g)\n" - " -m, --stress-ratio VAL applied stress ratio [-] (default %g)\n" - " -V, --velocity-bottom VAL base velocity at bottom [m/s] (default %g)\n" - " -A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default %g)\n" - " -b, --rate-dependence VAL rate dependence beyond yield [-] (default %g)\n" - " -f, --friction-coefficient VAL grain friction coefficient [-] (default %g)\n" - " -p, --porosity VAL porosity fraction [-] (default %g)\n" - " -d, --grain-size VAL representative grain size [m] (default %g)\n" - " -r, --density VAL grain material density [kg/m^3] (default %g)\n" - " -n, --resolution VAL number of cells in domain [-] (default %d)\n" - " -o, --origo VAL coordinate system origo [m] (default %g)\n" - " -L, --length VAL domain length [m] (default %g)\n" - " -F, --fluid enable pore fluid computations\n" - " -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default %g)\n" - " -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default %g)\n" - " -R, --fluid-density VAL fluid density [kg/m^3] (default %g)\n" - " -k, --fluid-permeability VAL fluid intrinsic permeability [m^2] (default %g)\n" - " -O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa] (default %g)\n" - " -a, --fluid-pressure-ampl VAL amplitude of pressure variations [Pa] (default %g)\n" - " -q, --fluid-pressure-freq VAL frequency of pressure variations [s^-1] (default %g)\n" - " -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa] (default %g)\n" - " -t, --time VAL simulation start time [s] (default %g)\n" - " -T, --time-end VAL simulation end time [s] (default %g)\n" - " -D, --time-step VAL computational time step length [s] (default %g)\n" - " -I, --file-interval VAL interval between output files [s] (default %g)\n" - " -v, --version show version information\n" - " -h, --help show this message\n", - __func__, PROGNAME, - sim.name, - sim.G, - sim.P_wall, - sim.mu_wall, - sim.v_x_bot, - sim.A, - sim.b, - sim.mu_s, - sim.phi[0], - sim.d, - sim.rho_s, - sim.nz, - sim.origo_z, - sim.L_z, - sim.beta_f, - sim.mu_f, - sim.rho_f, - sim.k[0], - sim.p_f_top, - sim.p_f_mod_ampl, - sim.p_f_mod_freq, - sim.p_f_mod_phase, - sim.t, - sim.t_end, - sim.dt, - sim.file_dt); - free(sim.phi); - free(sim.k); + struct simulation sim = init_sim(); + printf("%s: %s [OPTIONS] [NAME]\n" + "runs a simulation and outputs state in files prefixed with NAME.\n" + "If NAME is not specified, the default value '%s' is used.\n" + "optional arguments:\n" + " -N, --normalize normalize output velocity\n" + " -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n" + " -P, --normal-stress VAL normal stress on top [Pa] (default %g)\n" + " -m, --stress-ratio VAL applied stress ratio [-] (default %g)\n" + " -V, --velocity-bottom VAL base velocity at bottom [m/s] (default %g)\n" + " -A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default %g)\n" + " -b, --rate-dependence VAL rate dependence beyond yield [-] (default %g)\n" + " -f, --friction-coefficient VAL grain friction coefficient [-] (default %g)\n" + " -p, --porosity VAL porosity fraction [-] (default %g)\n" + " -d, --grain-size VAL representative grain size [m] (default %g)\n" + " -r, --density VAL grain material density [kg/m^3] (default %g)\n" + " -n, --resolution VAL number of cells in domain [-] (default %d)\n" + " -o, --origo VAL coordinate system origo [m] (default %g)\n" + " -L, --length VAL domain length [m] (default %g)\n" + " -F, --fluid enable pore fluid computations\n" + " -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default %g)\n" + " -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default %g)\n" + " -R, --fluid-density VAL fluid density [kg/m^3] (default %g)\n" + " -k, --fluid-permeability VAL fluid intrinsic permeability [m^2] (default %g)\n" + " -O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa] (default %g)\n" + " -a, --fluid-pressure-ampl VAL amplitude of pressure variations [Pa] (default %g)\n" + " -q, --fluid-pressure-freq VAL frequency of pressure variations [s^-1] (default %g)\n" + " -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa] (default %g)\n" + " -t, --time VAL simulation start time [s] (default %g)\n" + " -T, --time-end VAL simulation end time [s] (default %g)\n" + " -D, --time-step VAL computational time step length [s] (default %g)\n" + " -I, --file-interval VAL interval between output files [s] (default %g)\n" + " -v, --version show version information\n" + " -h, --help show this message\n", + __func__, PROGNAME, + sim.name, + sim.G, + sim.P_wall, + sim.mu_wall, + sim.v_x_bot, + sim.A, + sim.b, + sim.mu_s, + sim.phi[0], + sim.d, + sim.rho_s, + sim.nz, + sim.origo_z, + sim.L_z, + sim.beta_f, + sim.mu_f, + sim.rho_f, + sim.k[0], + sim.p_f_top, + sim.p_f_mod_ampl, + sim.p_f_mod_freq, + sim.p_f_mod_phase, + sim.t, + sim.t_end, + sim.dt, + sim.file_dt); + free(sim.phi); + free(sim.k); } static void version(void) { - printf("%s v%s\n" - "Licensed under the GNU Public License, v3+\n" - "written by Anders Damsgaard, anders@adamsgaard.dk\n" - "https://gitlab.com/admesg/1d_fd_simple_shear\n" - , PROGNAME, VERSION); + printf("%s v%s\n" + "Licensed under the GNU Public License, v3+\n" + "written by Anders Damsgaard, anders@adamsgaard.dk\n" + "https://gitlab.com/admesg/1d_fd_simple_shear\n" + , PROGNAME, VERSION); } int main(int argc, char* argv[]) { - /* load with default values */ - struct simulation sim = init_sim(); + /* load with default values */ + struct simulation sim = init_sim(); - int normalize = 0; + int normalize = 0; - int opt; - const char* optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T:D:I:"; - const struct option longopts[] = { - {"help", no_argument, NULL, 'h'}, - {"version", no_argument, NULL, 'v'}, - {"normalize", no_argument, NULL, 'N'}, - {"gravity", required_argument, NULL, 'G'}, - {"normal-stress", required_argument, NULL, 'P'}, - {"stress-ratio", required_argument, NULL, 'm'}, - {"velocity-bottom", required_argument, NULL, 'V'}, - {"nonlocal-amplitude", required_argument, NULL, 'A'}, - {"rate-dependence", required_argument, NULL, 'b'}, - {"friction-coefficient", required_argument, NULL, 'f'}, - {"porosity", required_argument, NULL, 'p'}, - {"grain-size", required_argument, NULL, 'd'}, - {"density", required_argument, NULL, 'r'}, - {"resolution", required_argument, NULL, 'n'}, - {"origo", required_argument, NULL, 'o'}, - {"length", required_argument, NULL, 'L'}, - {"fluid", no_argument, NULL, 'F'}, - {"fluid-compressiblity", required_argument, NULL, 'c'}, - {"fluid-viscosity", required_argument, NULL, 'i'}, - {"fluid-density", required_argument, NULL, 'R'}, - {"fluid-permeability", required_argument, NULL, 'k'}, - {"fluid-pressure-top", required_argument, NULL, 'O'}, - {"fluid-pressure-ampl", required_argument, NULL, 'a'}, - {"fluid-pressure-freq", required_argument, NULL, 'q'}, - {"fluid-pressure-phase", required_argument, NULL, 'H'}, - {"time", required_argument, NULL, 't'}, - {"time-end", required_argument, NULL, 'T'}, - {"time-step", required_argument, NULL, 'D'}, - {"file-interval", required_argument, NULL, 'I'}, - {NULL, 0, NULL, 0} - }; + int opt; + const char* optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T:D:I:"; + const struct option longopts[] = { + {"help", no_argument, NULL, 'h'}, + {"version", no_argument, NULL, 'v'}, + {"normalize", no_argument, NULL, 'N'}, + {"gravity", required_argument, NULL, 'G'}, + {"normal-stress", required_argument, NULL, 'P'}, + {"stress-ratio", required_argument, NULL, 'm'}, + {"velocity-bottom", required_argument, NULL, 'V'}, + {"nonlocal-amplitude", required_argument, NULL, 'A'}, + {"rate-dependence", required_argument, NULL, 'b'}, + {"friction-coefficient", required_argument, NULL, 'f'}, + {"porosity", required_argument, NULL, 'p'}, + {"grain-size", required_argument, NULL, 'd'}, + {"density", required_argument, NULL, 'r'}, + {"resolution", required_argument, NULL, 'n'}, + {"origo", required_argument, NULL, 'o'}, + {"length", required_argument, NULL, 'L'}, + {"fluid", no_argument, NULL, 'F'}, + {"fluid-compressiblity", required_argument, NULL, 'c'}, + {"fluid-viscosity", required_argument, NULL, 'i'}, + {"fluid-density", required_argument, NULL, 'R'}, + {"fluid-permeability", required_argument, NULL, 'k'}, + {"fluid-pressure-top", required_argument, NULL, 'O'}, + {"fluid-pressure-ampl", required_argument, NULL, 'a'}, + {"fluid-pressure-freq", required_argument, NULL, 'q'}, + {"fluid-pressure-phase", required_argument, NULL, 'H'}, + {"time", required_argument, NULL, 't'}, + {"time-end", required_argument, NULL, 'T'}, + {"time-step", required_argument, NULL, 'D'}, + {"file-interval", required_argument, NULL, 'I'}, + {NULL, 0, NULL, 0} + }; - double new_phi = sim.phi[0]; - double new_k = sim.k[0]; - while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) { - switch (opt) { - case -1: /* no more arguments */ - case 0: /* long options toggles */ - break; + double new_phi = sim.phi[0]; + double new_k = sim.k[0]; + while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) { + switch (opt) { + case -1: /* no more arguments */ + case 0: /* long options toggles */ + break; - case 'h': - usage(); - free(sim.phi); - free(sim.k); - return 0; - case 'v': - version(); - free(sim.phi); - free(sim.k); - return 0; - case 'n': - sim.nz = atoi(optarg); - break; - case 'N': - normalize = 1; - break; - case 'G': - sim.G = atof(optarg); - break; - case 'P': - sim.P_wall = atof(optarg); - break; - case 'm': - sim.mu_wall = atof(optarg); - break; - case 'V': - sim.v_x_bot = atof(optarg); - break; - case 'A': - sim.A = atof(optarg); - break; - case 'b': - sim.b = atof(optarg); - break; - case 'f': - sim.mu_s = atof(optarg); - break; - case 'p': - new_phi = atof(optarg); - break; - case 'd': - sim.d = atof(optarg); - break; - case 'r': - sim.rho_s = atof(optarg); - break; - case 'o': - sim.origo_z = atof(optarg); - break; - case 'L': - sim.L_z = atof(optarg); - break; - case 'F': - sim.fluid = 1; - break; - case 'c': - sim.beta_f = atof(optarg); - break; - case 'i': - sim.mu_f = atof(optarg); - break; - case 'R': - sim.rho_f = atof(optarg); - break; - case 'k': - new_k = atof(optarg); - break; - case 'O': - sim.p_f_top = atof(optarg); - break; - case 'a': - sim.p_f_mod_ampl = atof(optarg); - break; - case 'q': - sim.p_f_mod_freq = atof(optarg); - break; - case 'H': - sim.p_f_mod_phase = atof(optarg); - break; - case 't': - sim.t = atof(optarg); - break; - case 'T': - sim.t_end = atof(optarg); - break; - case 'D': - sim.dt = atof(optarg); - break; - case 'I': - sim.file_dt = atof(optarg); - break; + case 'h': + usage(); + free(sim.phi); + free(sim.k); + return 0; + case 'v': + version(); + free(sim.phi); + free(sim.k); + return 0; + case 'n': + sim.nz = atoi(optarg); + break; + case 'N': + normalize = 1; + break; + case 'G': + sim.G = atof(optarg); + break; + case 'P': + sim.P_wall = atof(optarg); + break; + case 'm': + sim.mu_wall = atof(optarg); + break; + case 'V': + sim.v_x_bot = atof(optarg); + break; + case 'A': + sim.A = atof(optarg); + break; + case 'b': + sim.b = atof(optarg); + break; + case 'f': + sim.mu_s = atof(optarg); + break; + case 'p': + new_phi = atof(optarg); + break; + case 'd': + sim.d = atof(optarg); + break; + case 'r': + sim.rho_s = atof(optarg); + break; + case 'o': + sim.origo_z = atof(optarg); + break; + case 'L': + sim.L_z = atof(optarg); + break; + case 'F': + sim.fluid = 1; + break; + case 'c': + sim.beta_f = atof(optarg); + break; + case 'i': + sim.mu_f = atof(optarg); + break; + case 'R': + sim.rho_f = atof(optarg); + break; + case 'k': + new_k = atof(optarg); + break; + case 'O': + sim.p_f_top = atof(optarg); + break; + case 'a': + sim.p_f_mod_ampl = atof(optarg); + break; + case 'q': + sim.p_f_mod_freq = atof(optarg); + break; + case 'H': + sim.p_f_mod_phase = atof(optarg); + break; + case 't': + sim.t = atof(optarg); + break; + case 'T': + sim.t_end = atof(optarg); + break; + case 'D': + sim.dt = atof(optarg); + break; + case 'I': + sim.file_dt = atof(optarg); + break; - default: - fprintf(stderr, "%s: invalid option -- %c\n", argv[0], opt); - fprintf(stderr, "Try `%s --help` for more information\n", - argv[0]); - return -2; - } - } - for (int i=optind; i<argc; ++i) { - if (i>optind) { - fprintf(stderr, "error: more than one simulation name specified\n"); - return 1; - } - sprintf(sim.name, "%s", argv[i]); - } + default: + fprintf(stderr, "%s: invalid option -- %c\n", argv[0], opt); + fprintf(stderr, "Try `%s --help` for more information\n", + argv[0]); + return -2; + } + } + for (int i=optind; i<argc; ++i) { + if (i>optind) { + fprintf(stderr, + "error: more than one simulation name specified\n"); + return 1; + } + sprintf(sim.name, "%s", argv[i]); + } - prepare_arrays(&sim); + prepare_arrays(&sim); - if (!isnan(new_phi)) - for (int i=0; i<sim.nz; ++i) - sim.phi[i] = new_phi; - if (!isnan(new_k)) - for (int i=0; i<sim.nz; ++i) - sim.k[i] = new_k; + if (!isnan(new_phi)) + for (int i=0; i<sim.nz; ++i) + sim.phi[i] = new_phi; + if (!isnan(new_k)) + for (int i=0; i<sim.nz; ++i) + sim.k[i] = new_k; - lithostatic_pressure_distribution(&sim); + lithostatic_pressure_distribution(&sim); - if (sim.fluid) - hydrostatic_fluid_pressure_distribution(&sim); + if (sim.fluid) + hydrostatic_fluid_pressure_distribution(&sim); #ifdef DEBUG - puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, sim.nz+2); - puts(""); - puts(".. normal stress before iterations:"); print_array(sim.sigma_n, sim.nz); - puts(""); + puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, sim.nz+2); + puts(""); + puts(".. normal stress before iterations:"); print_array(sim.sigma_n, sim.nz); + puts(""); #endif - double filetimeclock = 0.0; - unsigned long iter = 0; - while (sim.t <= sim.t_end) { + double filetimeclock = 0.0; + unsigned long iter = 0; + while (sim.t <= sim.t_end) { - if (sim.fluid) { - if (darcy_solver_1d(&sim, 10000, 1e-5)) - exit(1); + if (sim.fluid) { + if (darcy_solver_1d(&sim, 10000, 1e-5)) + exit(1); #ifdef DEBUG - puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2); - puts(""); + puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2); + puts(""); #endif - } + } - compute_effective_stress(&sim); - compute_friction(&sim); - compute_cooperativity_length(&sim); + compute_effective_stress(&sim); + compute_friction(&sim); + compute_cooperativity_length(&sim); - if (iter == 0) - check_simulation_parameters(&sim); + if (iter == 0) + check_simulation_parameters(&sim); #ifdef DEBUG - puts("\n## Before solver"); - puts(".. sigma_n_eff:"); print_array(sim.sigma_n_eff, sim.nz); - puts(".. mu:"); print_array(sim.mu, sim.nz); + puts("\n## Before solver"); + puts(".. sigma_n_eff:"); print_array(sim.sigma_n_eff, sim.nz); + puts(".. mu:"); print_array(sim.mu, sim.nz); #endif - if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5)) - exit(1); + if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5)) + exit(1); - compute_shear_strain_rate_plastic(&sim); - compute_shear_velocity(&sim); + compute_shear_strain_rate_plastic(&sim); + compute_shear_velocity(&sim); - sim.t += sim.dt; - filetimeclock += sim.dt; - iter++; + sim.t += sim.dt; + filetimeclock += sim.dt; + iter++; - if (filetimeclock >= sim.file_dt || iter == 0) { - write_output_file(&sim); - filetimeclock = 0.0; - } - } + if (filetimeclock >= sim.file_dt || iter == 0) { + write_output_file(&sim); + filetimeclock = 0.0; + } + } - if (normalize) { - double max_v_x = max(sim.v_x, sim.nz); - for (int i=0; i<sim.nz; ++i) - sim.v_x[i] /= max_v_x; - } + if (normalize) { + double max_v_x = max(sim.v_x, sim.nz); + for (int i=0; i<sim.nz; ++i) + sim.v_x[i] /= max_v_x; + } - if (sim.fluid) - for (int i=0; i<sim.nz; ++i) - printf("%.17g\t%.17g\t%.17g\t%.17g\n", - sim.z[i], - sim.v_x[i], - sim.sigma_n_eff[i], - sim.p_f_ghost[idx1g(i)]); - else - print_three_arrays(sim.z, sim.v_x, sim.sigma_n_eff, sim.nz); + if (sim.fluid) + for (int i=0; i<sim.nz; ++i) + printf("%.17g\t%.17g\t%.17g\t%.17g\n", + sim.z[i], + sim.v_x[i], + sim.sigma_n_eff[i], + sim.p_f_ghost[idx1g(i)]); + else + print_three_arrays(sim.z, sim.v_x, sim.sigma_n_eff, sim.nz); - free_arrays(&sim); - return 0; + free_arrays(&sim); + return 0; } (DIR) diff --git a/simulation.c b/simulation.c t@@ -6,430 +6,421 @@ void prepare_arrays(struct simulation* sim) { - 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->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 */ - free(sim->phi); - sim->phi = zeros(sim->nz); /* porosity */ - free(sim->k); - 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_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */ + 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->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 */ + free(sim->phi); + sim->phi = zeros(sim->nz); /* porosity */ + free(sim->k); + 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_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */ } void free_arrays(struct simulation* sim) { - free(sim->z); - free(sim->mu); - free(sim->sigma_n_eff); - free(sim->sigma_n); - free(sim->p_f_ghost); - free(sim->k); - free(sim->phi); - free(sim->xi); - free(sim->gamma_dot_p); - free(sim->v_x); - free(sim->g_ghost); + free(sim->z); + free(sim->mu); + free(sim->sigma_n_eff); + free(sim->sigma_n); + free(sim->p_f_ghost); + free(sim->k); + free(sim->phi); + free(sim->xi); + free(sim->gamma_dot_p); + free(sim->v_x); + free(sim->g_ghost); } static void warn_parameter_value( - const char message[], - const double value, - int* return_status) + const char message[], + const double value, + int* return_status) { - fprintf(stderr, "check_simulation_parameters: %s (%.17g)\n", - message, value); - *return_status = 1; + fprintf(stderr, "check_simulation_parameters: %s (%.17g)\n", + message, value); + *return_status = 1; } -static void check_float( - const char name[], - const double value, - int* return_status) +static void check_float(const char name[], + const double value, + int* return_status) { - if (isnan(value)) { - warn_parameter_value("%s is NaN", value, return_status); - *return_status = 1; - } else if (isinf(value)) { - warn_parameter_value("%s is infinite", value, return_status); - *return_status = 1; - } + if (isnan(value)) { + warn_parameter_value("%s is NaN", value, return_status); + *return_status = 1; + } else if (isinf(value)) { + warn_parameter_value("%s is infinite", value, return_status); + *return_status = 1; + } } void check_simulation_parameters(const struct simulation* sim) { - int return_status = 0; - - check_float("sim.G", sim->G, &return_status); - if (sim->G < 0.0) - warn_parameter_value("sim.G is negative", sim->G, &return_status); - - 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); - - 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); - - check_float("sim.A", sim->A, &return_status); - if (sim->A < 0.0) - warn_parameter_value("sim.A is negative", sim->A, &return_status); - - check_float("sim.b", sim->b, &return_status); - if (sim->b < 0.0) - warn_parameter_value("sim.b is negative", sim->b, &return_status); - - 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); - - 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); - - 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); - - if (sim->nz <= 0) - warn_parameter_value("sim.nz is not a positive number", sim->nz, - &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); - - if (sim->nz <= 0) - warn_parameter_value("sim.nz is not a positive number", sim->nz, - &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); - - 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); - - 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); - - check_float("sim.dt", sim->t_end, &return_status); - if (sim->dt <= 0.0) - warn_parameter_value("sim.dt is not a positive number", - 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); - - 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); - - if (sim->fluid != 0 && sim->fluid != 1) - warn_parameter_value("sim.fluid has an invalid value", - (double)sim->fluid, &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); - - 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); - - 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); - - 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); - - 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); - - 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 (return_status != 0) { - fprintf(stderr, "error: aborting due to invalid parameter choices\n"); - exit(return_status); - } + int return_status = 0; + + check_float("sim.G", sim->G, &return_status); + if (sim->G < 0.0) + warn_parameter_value("sim.G is negative", sim->G, &return_status); + + 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); + + 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); + + check_float("sim.A", sim->A, &return_status); + if (sim->A < 0.0) + warn_parameter_value("sim.A is negative", sim->A, &return_status); + + check_float("sim.b", sim->b, &return_status); + if (sim->b < 0.0) + warn_parameter_value("sim.b is negative", sim->b, &return_status); + + 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); + + 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); + + 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); + + if (sim->nz <= 0) + warn_parameter_value("sim.nz is not a positive number", sim->nz, + &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); + + if (sim->nz <= 0) + warn_parameter_value("sim.nz is not a positive number", sim->nz, + &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); + + 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); + + 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); + + check_float("sim.dt", sim->t_end, &return_status); + if (sim->dt <= 0.0) + warn_parameter_value("sim.dt is not a positive number", + 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); + + 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); + + if (sim->fluid != 0 && sim->fluid != 1) + warn_parameter_value("sim.fluid has an invalid value", + (double)sim->fluid, &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); + + 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); + + 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); + + 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); + + 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); + + 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 (return_status != 0) { + fprintf(stderr, "error: aborting due to invalid parameter choices\n"); + exit(return_status); + } } void lithostatic_pressure_distribution(struct simulation* sim) { - for (int 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]); + for (int 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]); } void compute_friction(struct simulation* sim) { - if (sim->fluid) - for (int 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 (int 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); - + if (sim->fluid) + for (int 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 (int 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); } double shear_strain_rate_plastic( - const double fluidity, - const double friction) + const double fluidity, + const double friction) { - return fluidity*friction; + return fluidity*friction; } void compute_shear_strain_rate_plastic(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) - sim->gamma_dot_p[i] = - shear_strain_rate_plastic(sim->g_ghost[idx1g(i)], sim->mu[i]); + for (int i=0; i<sim->nz; ++i) + sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[idx1g(i)], + sim->mu[i]); } void compute_shear_velocity(struct simulation* sim) { - // TODO: implement iterative solver - // Dirichlet BC at bottom - sim->v_x[0] = sim->v_x_bot; + // TODO: implement iterative solver + // Dirichlet BC at bottom + sim->v_x[0] = sim->v_x_bot; - for (int i=1; i<sim->nz; ++i) - sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz; + for (int 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) { - if (sim->fluid) - for (int i=0; i<sim->nz; ++i) - sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)]; - else - for (int i=0; i<sim->nz; ++i) - sim->sigma_n_eff[i] = sim->sigma_n[i]; + if (sim->fluid) + for (int i=0; i<sim->nz; ++i) + sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)]; + else + for (int i=0; i<sim->nz; ++i) + sim->sigma_n_eff[i] = sim->sigma_n[i]; } -double cooperativity_length( - const double A, - const double d, - const double mu, - const double mu_s) +double cooperativity_length(const double A, + const double d, + const double mu, + const double mu_s) { - return A*d/sqrt(fabs(mu - mu_s)); + return A*d/sqrt(fabs(mu - mu_s)); } void compute_cooperativity_length(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) - sim->xi[i] = cooperativity_length( - sim->A, sim->d, sim->mu[i], sim->mu_s); + for (int i=0; i<sim->nz; ++i) + sim->xi[i] = cooperativity_length(sim->A, sim->d, sim->mu[i], + sim->mu_s); } -double local_fluidity( - const double p, - const double mu, - const double mu_s, - const double b, - const double rho_s, - const double d) +double local_fluidity(const double p, + const double mu, + const double mu_s, + const double b, + const double rho_s, + const double d) { - if (mu <= mu_s) - return 0.0; - else - return sqrt(p/rho_s*d*d) * (mu - mu_s)/(b*mu); + if (mu <= mu_s) + return 0.0; + else + return sqrt(p/rho_s*d*d) * (mu - mu_s)/(b*mu); } void compute_local_fluidity(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) - sim->g_ghost[idx1g(i)] = local_fluidity( - sim->sigma_n_eff[i], - sim->mu[i], - sim->mu_s, - sim->b, - sim->rho_s, - sim->d); + for (int i=0; i<sim->nz; ++i) + sim->g_ghost[idx1g(i)] = local_fluidity(sim->sigma_n_eff[i], + sim->mu[i], + sim->mu_s, + sim->b, + sim->rho_s, + sim->d); } void set_bc_neumann(double* g_ghost, const int nz, const int boundary) { - if (boundary == -1) - g_ghost[0] = g_ghost[1]; - else if (boundary == +1) - g_ghost[nz+1] = g_ghost[nz]; - else { - fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary); - exit(1); - } + if (boundary == -1) + g_ghost[0] = g_ghost[1]; + else if (boundary == +1) + g_ghost[nz+1] = g_ghost[nz]; + else { + fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary); + exit(1); + } } -void set_bc_dirichlet( - double* g_ghost, - const int nz, - const int boundary, - const double value) +void set_bc_dirichlet(double* g_ghost, + const int nz, + const int boundary, + const double value) { - if (boundary == -1) - g_ghost[0] = value; - else if (boundary == +1) - g_ghost[nz+1] = value; - else { - fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary); - exit(1); - } + if (boundary == -1) + g_ghost[0] = value; + else if (boundary == +1) + g_ghost[nz+1] = value; + else { + fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary); + exit(1); + } } -void poisson_solver_1d_cell_update( - int i, - const double* g_in, - double* g_out, - double* r_norm, - const double dz, - const double* mu, - const double* p, - const double* xi, - const double mu_s, - const double b, - const double rho_s, - const double d) +void poisson_solver_1d_cell_update(int i, + const double* g_in, + double* g_out, + double* r_norm, + const double dz, + const double* mu, + const double* p, + const double* xi, + const double mu_s, + const double b, + const double rho_s, + const double d) { - double coorp_term = dz*dz/(2.0*pow(xi[i], 2.0)); - int gi = idx1g(i); - g_out[gi] = 1.0/(1.0 + coorp_term)*(coorp_term* - local_fluidity(p[i], mu[i], mu_s, b, rho_s, d) - + g_in[gi+1]/2.0 - + g_in[gi-1]/2.0); + double coorp_term = dz*dz/(2.0*pow(xi[i], 2.0)); + int gi = idx1g(i); + g_out[gi] = 1.0/(1.0 + coorp_term)*(coorp_term* + local_fluidity(p[i], mu[i], mu_s, b, rho_s, d) + + g_in[gi+1]/2.0 + + g_in[gi-1]/2.0); - r_norm[i] = pow(g_out[gi] - g_in[gi], 2.0) / (pow(g_out[gi], 2.0) + 1e-16); + r_norm[i] = pow(g_out[gi] - g_in[gi], 2.0) / (pow(g_out[gi], 2.0) + 1e-16); #ifdef DEBUG - printf("-- %d --------------\n", i); - printf("coorp_term: %g\n", coorp_term); - printf(" g_local: %g\n", local_fluidity(p[i], mu[i], mu_s, b, rho_s, d)); - printf(" g_in: %g\n", g_in[gi]); - printf(" g_out: %g\n", g_out[gi]); - printf(" r_norm: %g\n", r_norm[i]); + printf("-- %d --------------\n", i); + printf("coorp_term: %g\n", coorp_term); + printf(" g_local: %g\n", local_fluidity(p[i], mu[i], mu_s, b, rho_s, d)); + printf(" g_in: %g\n", g_in[gi]); + printf(" g_out: %g\n", g_out[gi]); + printf(" r_norm: %g\n", r_norm[i]); #endif } -int implicit_1d_jacobian_poisson_solver( - struct simulation* sim, - const int max_iter, - const double rel_tol) +int implicit_1d_jacobian_poisson_solver(struct simulation* sim, + const int max_iter, + const double rel_tol) { - double* g_ghost_out = empty(sim->nz+2); - double* r_norm = empty(sim->nz); + double* g_ghost_out = empty(sim->nz+2); + double* r_norm = empty(sim->nz); - int iter; - double r_norm_max = NAN; + int iter; + double r_norm_max = NAN; - for (iter=0; iter<max_iter; ++iter) { + for (iter=0; iter<max_iter; ++iter) { #ifdef DEBUG - printf("\n@@@ ITERATION %d @@@\n", iter); + printf("\n@@@ ITERATION %d @@@\n", iter); #endif - /* Dirichlet BCs resemble fixed particle velocities */ - set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); - set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0); - - /* Neumann BCs resemble free surfaces */ - /* set_bc_neumann(sim->g_ghost, sim->nz, +1); */ - - for (int i=0; i<sim->nz; ++i) - poisson_solver_1d_cell_update( - i, - sim->g_ghost, - g_ghost_out, - r_norm, - sim->dz, - sim->mu, - sim->sigma_n_eff, - sim->xi, - sim->mu_s, - sim->b, - sim->rho_s, - sim->d); - r_norm_max = max(r_norm, sim->nz); - - 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); */ - return 0; - } - } - - free(g_ghost_out); - free(r_norm); - fprintf(stderr, "implicit_1d_jacobian_poisson_solver: "); - fprintf(stderr, "Solution did not converge after %d iterations\n", iter); - fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); - return 1; + /* Dirichlet BCs resemble fixed particle velocities */ + set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0); + set_bc_dirichlet(sim->g_ghost, sim->nz, +1, 0.0); + + /* Neumann BCs resemble free surfaces */ + /* set_bc_neumann(sim->g_ghost, sim->nz, +1); */ + + for (int i=0; i<sim->nz; ++i) + poisson_solver_1d_cell_update(i, + sim->g_ghost, + g_ghost_out, + r_norm, + sim->dz, + sim->mu, + sim->sigma_n_eff, + sim->xi, + sim->mu_s, + sim->b, + sim->rho_s, + sim->d); + r_norm_max = max(r_norm, sim->nz); + + 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); */ + return 0; + } + } + + free(g_ghost_out); + free(r_norm); + fprintf(stderr, "implicit_1d_jacobian_poisson_solver: "); + fprintf(stderr, "Solution did not converge after %d iterations\n", iter); + fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max); + return 1; } void write_output_file(struct simulation* sim) { - char outfile[200]; - FILE *fp; - sprintf(outfile, "%s.output%05d.txt", sim->name, sim->n_file++); - - fp = fopen(outfile, "w"); - if (sim->fluid) - for (int i=0; i<sim->nz; ++i) - fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n", - sim->z[i], - sim->v_x[i], - sim->sigma_n_eff[i], - sim->p_f_ghost[idx1g(i)], - sim->mu[i]); - else - fprint_three_arrays(fp, sim->z, sim->v_x, sim->sigma_n_eff, sim->nz); - - fclose(fp); + char outfile[200]; + FILE *fp; + sprintf(outfile, "%s.output%05d.txt", sim->name, sim->n_file++); + + fp = fopen(outfile, "w"); + if (sim->fluid) + for (int i=0; i<sim->nz; ++i) + fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n", + sim->z[i], + sim->v_x[i], + sim->sigma_n_eff[i], + sim->p_f_ghost[idx1g(i)], + sim->mu[i]); + else + fprint_three_arrays(fp, sim->z, sim->v_x, sim->sigma_n_eff, sim->nz); + + fclose(fp); }