cngf-pf

continuum model for granular flows with pore-pressure dynamics
git clone git://src.adamsgaard.dk/1d_fd_simple_shear
Log | Files | Refs | README | LICENSE Back to index

1d_fd_simple_shear.c (5723B)


      1 #include <stdlib.h>
      2 #include <math.h>
      3 #include <string.h>
      4 #include <time.h>
      5 #include <unistd.h>
      6 
      7 #include "simulation.h"
      8 #include "fluid.h"
      9 
     10 #include "arg.h"
     11 
     12 /* relative tolerance criteria for granular fluidity solver */
     13 #define RTOL 1e-5
     14 #define MAX_ITER_1D_FD_SIMPLE_SHEAR 100000
     15 
     16 /* uncomment to print time spent per time step to stdout */
     17 /* #define BENCHMARK_PERFORMANCE */
     18 
     19 char *argv0;
     20 
     21 static void
     22 usage(void)
     23 {
     24 	fprintf(stderr, "usage: %s "
     25 		"[-A grain-nonlocal-ampl] "
     26 		"[-a fluid-pressure-ampl] "
     27 		"[-b grain-rate-dependence] "
     28 		"[-C fluid-compressibility] "
     29 		"[-c grain-cohesion] "
     30 		"[-d grain-size] "
     31 		"[-D fluid-diffusivity] "
     32 		"[-e end-time] "
     33 		"[-F] "
     34 		"[-f applied-shear-friction] "
     35 		"[-g gravity-accel] "
     36 		"[-H fluid-pressure-phase] "
     37 		"[-h] "
     38 		"[-I file-interval] "
     39 		"[-i fluid-viscosity] "
     40 		"[-K dilatancy-constant] "
     41 		"[-k fluid-permeability] "
     42 		"[-L length] "
     43 		"[-l applied-shear-vel-limit] "
     44 		"[-m grain-friction] "
     45 		"[-N] "
     46 		"[-n normal-stress] "
     47 		"[-O fluid-pressure-top] "
     48 		"[-o origo] "
     49 		"[-P grain-compressibility] "
     50 		"[-p grain-porosity] "
     51 		"[-q fluid-pressure-freq] "
     52 		"[-R fluid-density] "
     53 		"[-r grain-density] "
     54 		"[-S fluid-pressure-pulse-shape] "
     55 		"[-s applied-shear-vel] "
     56 		"[-T] "
     57 		"[-t curr-time] "
     58 		"[-U resolution] "
     59 		"[-u fluid-pulse-time] "
     60 		"[-v] "
     61 		"[-Y max-porosity] "
     62 		"[-y min-porosity] "
     63 		"[name]\n", argv0);
     64 	exit(1);
     65 }
     66 
     67 int
     68 main(int argc, char *argv[])
     69 {
     70 	int i, normalize;
     71 	unsigned long iter;
     72 	double new_phi, new_k, filetimeclock;
     73 	struct simulation sim;
     74 
     75 #ifdef BENCHMARK_PERFORMANCE
     76 	clock_t t_begin, t_end;
     77 	double t_elapsed;
     78 
     79 #endif
     80 
     81 #ifdef __OpenBSD__
     82 	if (pledge("stdio wpath cpath", NULL) == -1) {
     83 		fprintf(stderr, "error: pledge failed");
     84 		exit(2);
     85 	}
     86 #endif
     87 
     88 	init_sim(&sim);
     89 	normalize = 0;
     90 	new_phi = sim.phi[0];
     91 	new_k = sim.k[0];
     92 
     93 	ARGBEGIN {
     94 	case 'A':
     95 		sim.A = atof(EARGF(usage()));
     96 		break;
     97 	case 'a':
     98 		sim.p_f_mod_ampl = atof(EARGF(usage()));
     99 		break;
    100 	case 'b':
    101 		sim.b = atof(EARGF(usage()));
    102 		break;
    103 	case 'C':
    104 		sim.beta_f = atof(EARGF(usage()));
    105 		break;
    106 	case 'c':
    107 		sim.C = atof(EARGF(usage()));
    108 		break;
    109 	case 'd':
    110 		sim.d = atof(EARGF(usage()));
    111 		break;
    112 	case 'D':
    113 		sim.D = atof(EARGF(usage()));
    114 		break;
    115 	case 'e':
    116 		sim.t_end = atof(EARGF(usage()));
    117 		break;
    118 	case 'F':
    119 		sim.fluid = 1;
    120 		break;
    121 	case 'f':
    122 		sim.mu_wall = atof(EARGF(usage()));
    123 		break;
    124 	case 'g':
    125 		sim.G = atof(EARGF(usage()));
    126 		break;
    127 	case 'H':
    128 		sim.p_f_mod_phase = atof(EARGF(usage()));
    129 		break;
    130 	case 'h':
    131 		usage();
    132 		break;
    133 	case 'I':
    134 		sim.file_dt = atof(EARGF(usage()));
    135 		break;
    136 	case 'i':
    137 		sim.mu_f = atof(EARGF(usage()));
    138 		break;
    139 	case 'K':
    140 		sim.dilatancy_angle = atof(EARGF(usage()));
    141 		break;
    142 	case 'k':
    143 		new_k = atof(EARGF(usage()));
    144 		break;
    145 	case 'L':
    146 		sim.L_z = atof(EARGF(usage()));
    147 		break;
    148 	case 'l':
    149 		sim.v_x_limit = atof(EARGF(usage()));
    150 		break;
    151 	case 'm':
    152 		sim.mu_s = atof(EARGF(usage()));
    153 		break;
    154 	case 'N':
    155 		normalize = 1;
    156 		break;
    157 	case 'n':
    158 		sim.P_wall = atof(EARGF(usage()));
    159 		break;
    160 	case 'O':
    161 		sim.p_f_top = atof(EARGF(usage()));
    162 		break;
    163 	case 'o':
    164 		sim.origo_z = atof(EARGF(usage()));
    165 		break;
    166 	case 'P':
    167 		sim.alpha = atof(EARGF(usage()));
    168 		break;
    169 	case 'p':
    170 		new_phi = atof(EARGF(usage()));
    171 		break;
    172 	case 'q':
    173 		sim.p_f_mod_freq = atof(EARGF(usage()));
    174 		break;
    175 	case 'R':
    176 		sim.rho_f = atof(EARGF(usage()));
    177 		break;
    178 	case 'r':
    179 		sim.rho_s = atof(EARGF(usage()));
    180 		break;
    181 	case 'S':
    182 		if (argv[1] == NULL)
    183 			usage();
    184 		else if (!strncmp(argv[1], "triangle", 8))
    185 			sim.p_f_mod_pulse_shape = 0;
    186 		else if (!strncmp(argv[1], "square", 6))
    187 			sim.p_f_mod_pulse_shape = 1;
    188 		else {
    189 			fprintf(stderr, "error: invalid pulse shape '%s'\n",
    190 				argv[1]);
    191 			return 1;
    192 		}
    193 		argc--;
    194 		argv++;
    195 		break;
    196 	case 's':
    197 		sim.v_x_fix = atof(EARGF(usage()));
    198 		break;
    199 	case 'T':
    200 		sim.transient = 1;
    201 		break;
    202 	case 't':
    203 		sim.t = atof(EARGF(usage()));
    204 		break;
    205 	case 'U':
    206 		sim.nz = atoi(EARGF(usage()));
    207 		break;
    208 	case 'u':
    209 		sim.p_f_mod_pulse_time = atof(EARGF(usage()));
    210 		break;
    211 	case 'V':
    212 		sim.v_x_bot = atof(EARGF(usage()));
    213 		break;
    214 	case 'v':
    215 		printf("%s-" VERSION "\n", argv0);
    216 		return 0;
    217 	case 'Y':
    218 		sim.phi_max = atof(EARGF(usage()));
    219 		break;
    220 	case 'y':
    221 		sim.phi_min = atof(EARGF(usage()));
    222 		break;
    223 	default:
    224 		usage();
    225 	} ARGEND;
    226 
    227 	if (argc == 1 && argv[0])
    228 		snprintf(sim.name, sizeof(sim.name), "%s", argv[0]);
    229 	else if (argc > 1)
    230 		usage();
    231 
    232 	if (sim.nz < 1)
    233 		sim.nz = (int) ceil(sim.L_z / sim.d);
    234 
    235 	prepare_arrays(&sim);
    236 
    237 	if (!isnan(new_phi))
    238 		for (i = 0; i < sim.nz; ++i)
    239 			sim.phi[i] = new_phi;
    240 	if (!isnan(new_k))
    241 		for (i = 0; i < sim.nz; ++i)
    242 			sim.k[i] = new_k;
    243 
    244 	lithostatic_pressure_distribution(&sim);
    245 
    246 	if (sim.fluid) {
    247 		hydrostatic_fluid_pressure_distribution(&sim);
    248 		if (set_largest_fluid_timestep(&sim, 0.5)) {
    249 			free_arrays(&sim);
    250 			return 20;
    251 		}
    252 	}
    253 	if (sim.dt > sim.file_dt)
    254 		sim.dt = sim.file_dt;
    255 	compute_effective_stress(&sim);
    256 
    257 	check_simulation_parameters(&sim);
    258 
    259 	filetimeclock = 0.0;
    260 	iter = 0;
    261 	do {
    262 
    263 #ifdef BENCHMARK_PERFORMANCE
    264 		t_begin = clock();
    265 #endif
    266 
    267 		if (coupled_shear_solver(&sim, MAX_ITER_1D_FD_SIMPLE_SHEAR, RTOL)) {
    268 			free_arrays(&sim);
    269 			exit(10);
    270 		}
    271 #ifdef BENCHMARK_PERFORMANCE
    272 		t_end = clock();
    273 		t_elapsed = (double) (t_end - t_begin) / CLOCKS_PER_SEC;
    274 		printf("time spent per time step = %g s\n", t_elapsed);
    275 #endif
    276 
    277 		if ((filetimeclock >= sim.file_dt || iter == 1) &&
    278 		    strncmp(sim.name, DEFAULT_SIMULATION_NAME,
    279 			    sizeof(DEFAULT_SIMULATION_NAME)) != 0) {
    280 			write_output_file(&sim, normalize);
    281 			filetimeclock = 0.0;
    282 		}
    283 		sim.t += sim.dt;
    284 		filetimeclock += sim.dt;
    285 		iter++;
    286 
    287 	} while (sim.t - sim.dt < sim.t_end);
    288 
    289 	print_output(&sim, stdout, normalize);
    290 
    291 	free_arrays(&sim);
    292 	return 0;
    293 }