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