simulation.h (6534B)
1 #ifndef SIMULATION_ 2 #define SIMULATION_ 3 4 #include "arrays.h" 5 6 #define DEFAULT_SIMULATION_NAME "unnamed_simulation" 7 #define PI 3.14159265358979323846 8 #define DEG2RAD(x) (x * PI / 180.0) 9 10 #ifdef __cplusplus 11 extern "C" { 12 #endif 13 14 /* Simulation settings */ 15 struct simulation { 16 17 /* simulation name to use for output files */ 18 char name[100]; 19 20 /* gravitational acceleration magnitude [m/s^2] */ 21 double G; 22 23 /* normal stress from the top wall [Pa] */ 24 double P_wall; 25 26 /* optionally fix top shear velocity to this value [m/s] */ 27 double v_x_fix; 28 29 /* optionally fix top shear velocity to this value [m/s] */ 30 double v_x_limit; 31 32 /* bottom velocity along x [m/s] */ 33 double v_x_bot; 34 35 /* stress ratio at top wall */ 36 double mu_wall; 37 38 /* nonlocal amplitude [-] */ 39 double A; 40 41 /* rate dependence beyond yield [-] */ 42 double b; 43 44 /* bulk and critical state static yield friction coefficient [-] */ 45 double mu_s; 46 47 /* material cohesion [Pa] */ 48 double C; 49 50 /* representative grain size [m] */ 51 double d; /* ohlala */ 52 53 /* grain material density [kg/m^3] */ 54 double rho_s; 55 56 /* nodes along z */ 57 int nz; 58 59 /* origo of axis [m] */ 60 double origo_z; 61 62 /* length of domain [m] */ 63 double L_z; 64 65 /* array of cell coordinates */ 66 double *z; 67 68 /* cell spacing [m] */ 69 double dz; 70 71 /* current time [s] */ 72 double t; 73 74 /* end time [s] */ 75 double t_end; 76 77 /* time step length [s] */ 78 double dt; 79 80 /* interval between output files [s] */ 81 double file_dt; 82 83 /* output file number */ 84 int n_file; 85 86 double transient; 87 double phi_min; 88 double phi_max; 89 double dilatancy_constant; 90 91 /* Fluid parameters */ 92 int fluid; /* flag to switch fluid on (1) or off (0) */ 93 double p_f_top; /* fluid pressure at the top [Pa] */ 94 double p_f_mod_ampl; /* amplitude of fluid pressure variations [Pa] */ 95 double p_f_mod_freq; /* frequency of fluid pressure variations [s^-1] */ 96 double p_f_mod_phase; /* phase of fluid pressure variations [s^-1] */ 97 double p_f_mod_pulse_time; /* single pressure pulse at this time [s] */ 98 int p_f_mod_pulse_shape; /* waveform for fluid-pressure pulse */ 99 double beta_f; /* adiabatic fluid compressibility [Pa^-1] */ 100 double alpha; /* adiabatic grain compressibility [Pa^-1] */ 101 double mu_f; /* fluid dynamic viscosity [Pa*s] */ 102 double rho_f; /* fluid density [kg/m^3] */ 103 double D; /* diffusivity [m^2/s], overrides k, beta_f, alpha, mu_f */ 104 105 /* arrays */ 106 double *mu; /* static yield friction [-] */ 107 double *mu_c; /* critical-state static yield friction [-] */ 108 double *sigma_n_eff; /* effective normal pressure [Pa] */ 109 double *sigma_n; /* normal stress [Pa] */ 110 double *p_f_ghost; /* fluid pressure [Pa] */ 111 double *p_f_next_ghost; /* fluid pressure for next iteration [Pa] */ 112 double *p_f_dot; /* fluid pressure change [Pa/s] */ 113 double *p_f_dot_impl; /* fluid pressure change (implicit solution) [Pa/s] */ 114 115 double *k; /* hydraulic permeability [m^2] */ 116 double *phi; /* porosity [-] */ 117 double *phi_c; /* critical-state porosity [-] */ 118 double *phi_dot; /* porosity change [s^-1] */ 119 double *xi; /* cooperativity length */ 120 double *gamma_dot_p; /* plastic shear strain rate [s^-1] */ 121 double *v_x; /* shear velocity [m/s] */ 122 double *d_x; /* cumulative shear displacement [m] */ 123 double *g_local; /* local fluidity */ 124 double *g_ghost; /* fluidity with ghost nodes */ 125 double *g_r_norm; /* normalized residual of fluidity field */ 126 double *I; /* inertia number [-] */ 127 double *tan_psi; /* tan(dilatancy_angle) [-] */ 128 double *old_val; /* temporary storage for iterative solvers */ 129 130 /* Persistent solver workspace (size nz unless noted otherwise), allocated 131 * once in prepare_arrays() and reused by Darcy/Poisson TDMA paths. */ 132 double *tdma_a; /* shared TDMA sub-diagonal coefficients */ 133 double *tdma_b; /* shared TDMA diagonal coefficients */ 134 double *tdma_c; /* shared TDMA super-diagonal coefficients */ 135 double *tdma_d; /* shared TDMA RHS coefficients */ 136 double *tdma_x; /* shared TDMA solution vector */ 137 double *tdma_c_prime; /* shared TDMA forward-sweep scratch */ 138 double *tdma_d_prime; /* shared TDMA forward-sweep scratch */ 139 double *darcy_k_n; /* Darcy predictor permeability workspace */ 140 double *darcy_phi_n; /* Darcy predictor porosity workspace */ 141 }; 142 143 void init_sim(struct simulation *sim); 144 145 /* returns: 0 ok, 1 grid size (nz) is less than 2 */ 146 int prepare_arrays(struct simulation *sim); 147 148 void free_arrays(struct simulation *sim); 149 150 /* Reset per-column state (stresses, velocities, fluidity, time, output-file 151 * counter) without reallocating; geometry (nz, L_z, z, dz), parameters, and 152 * solver workspaces are kept. phi and k are restored to the init_sim() 153 * defaults; set per-column values after the reset. For reuse of one prepared 154 * simulation across many independent columns. */ 155 void reset_column(struct simulation *sim); 156 157 /* returns: 0 ok, 1 one or more parameters are invalid */ 158 int check_simulation_parameters(struct simulation *sim); 159 void lithostatic_pressure_distribution(struct simulation *sim); 160 void compute_effective_stress(struct simulation *sim); 161 162 void set_bc_neumann(double *a, const int nz, const int boundary, 163 const double df, const double dx); 164 165 void set_bc_dirichlet(double *a, const int nz, const int boundary, 166 const double value); 167 168 double residual(double new_val, double old_val); 169 double kozeny_carman(const double diameter, const double porosity); 170 171 void write_output_file(struct simulation *sim, const int normalize); 172 void print_output(struct simulation *sim, FILE *fp, const int normalize); 173 174 /* returns: 0 ok, 1 transient solution not converged, 10 stress solution not 175 * converged, 11 fluid (Darcy) solver failed, 12 fluidity (Poisson) solver 176 * failed, 13 NaN in temporal increment */ 177 int coupled_shear_solver(struct simulation *sim, const int max_iter, 178 const double rel_tol); 179 180 void set_coupled_fluid_transient_timestep(struct simulation *sim, 181 const double safety); 182 183 double find_flux(const struct simulation *sim); 184 185 void print_solver_stats(FILE *fp); 186 void reset_solver_stats(void); 187 188 void tridiagonal_solver(double *x, const double *a, const double *b, 189 const double *c, const double *d, double *c_prime, 190 double *d_prime, int n); 191 192 #ifdef __cplusplus 193 } 194 #endif 195 196 #endif