simulation.h (5250B)
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 double *p_f_dot_impl_r_norm; /* normalized residual fluid pressure change [-] */ 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 double *fluid_old_val;/* temporary storage for fluid iterative solver */ 129 double *tmp_ghost; /* temporary storage for iterative solvers */ 130 double *p_f_dot_old; /* temporary storage for old p_f_dot */ 131 }; 132 133 void init_sim(struct simulation *sim); 134 void prepare_arrays(struct simulation *sim); 135 void free_arrays(struct simulation *sim); 136 void check_simulation_parameters(struct simulation *sim); 137 void lithostatic_pressure_distribution(struct simulation *sim); 138 void compute_effective_stress(struct simulation *sim); 139 140 void set_bc_neumann(double *a, 141 const int nz, 142 const int boundary, 143 const double df, 144 const double dx); 145 146 void set_bc_dirichlet(double *a, 147 const int nz, 148 const int boundary, 149 const double value); 150 151 double residual(double new_val, double old_val); 152 double kozeny_carman(const double diameter, const double porosity); 153 154 void write_output_file(struct simulation *sim, const int normalize); 155 void print_output(struct simulation *sim, FILE *fp, const int normalize); 156 157 int coupled_shear_solver(struct simulation *sim, 158 const int max_iter, 159 const double rel_tol); 160 161 void set_coupled_fluid_transient_timestep(struct simulation *sim, const double safety); 162 163 double find_flux(const struct simulation *sim); 164 165 #endif