tAdd makefile target for cyclic fluid test, fix k values at boundary, switchabl… | |
git clone git://src.adamsgaard.dk/cngf-pf | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
commit 1157d4b10b28dcb5291206da4dbab7a410f1fd10 | |
parent fc8db4ee3cef05e5f3620fc0f62c6733fb61ad38 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 16 Apr 2019 17:09:06 +0200 | |
Add makefile target for cyclic fluid test, fix k values at boundary, switchable… | |
Diffstat: | |
M 1d_fd_simple_shear.png | 0 | |
A 1d_fd_simple_shear_fluid.gp | 25 +++++++++++++++++++++++++ | |
M 1d_fd_simple_shear_rheology.png | 0 | |
M 1d_fd_simple_shear_rheology_iverso… | 0 | |
M 1d_fd_simple_shear_rheology_kamb.p… | 0 | |
M 1d_fd_simple_shear_rheology_tulacz… | 0 | |
M Makefile | 9 +++++++++ | |
M fluid.c | 57 ++++++++++++++++++++++-------… | |
M main.c | 21 +++++++++++++-------- | |
M simulation.c | 3 ++- | |
10 files changed, 89 insertions(+), 26 deletions(-) | |
--- | |
diff --git a/1d_fd_simple_shear.png b/1d_fd_simple_shear.png | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_fluid.gp b/1d_fd_simple_shear_fluid.gp | |
t@@ -0,0 +1,25 @@ | |
+#!/usr/bin/env gnuplot | |
+ | |
+# specify input file with: | |
+# gnuplot -e "filename='file.txt'" 1d_fd_simple_shear_fluid.gp | |
+ | |
+set terminal pngcairo color size 60 cm, 17.6 cm | |
+set multiplot layout 1,3 | |
+ | |
+set ylabel "Vertical position, z [m]" offset 2 | |
+ | |
+set yrange [0.0:0.73] | |
+#set xrange [-0.1:1.05] | |
+ | |
+set key bottom right #samplen 0.9 | |
+ | |
+set xlabel "Normalized horizontal velocity, v_x [-]" | |
+plot filename u 2:1 w l lw 2 lc "red" t "" | |
+ | |
+set xlabel "Effective normal stress [kPa]" | |
+plot filename u ($3/1000):1 w l lw 2 lc "black" t "" | |
+ | |
+set xlabel "Water pressure, p_f [kPa]" | |
+plot filename u ($4/1000):1 w l lw 2 lc "blue" t "" | |
+ | |
+unset multiplot | |
diff --git a/1d_fd_simple_shear_rheology.png b/1d_fd_simple_shear_rheology.png | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_rheology_iverson.png b/1d_fd_simple_shear_rheol… | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_rheology_kamb.png b/1d_fd_simple_shear_rheology… | |
Binary files differ. | |
diff --git a/1d_fd_simple_shear_rheology_tulaczyk.png b/1d_fd_simple_shear_rheo… | |
Binary files differ. | |
diff --git a/Makefile b/Makefile | |
t@@ -14,6 +14,15 @@ plots: 1d_fd_simple_shear.png \ | |
1d_fd_simple_shear_rheology_iverson.png \ | |
1d_fd_simple_shear_rheology_tulaczyk.png | |
+.PHONY: diurnal | |
+diurnaltest: 1d_fd_simple_shear 1d_fd_simple_shear_fluid.gp | |
+ /usr/bin/env zsh -c '\ | |
+ ./$< --resolution 10 --normal-stress 150e3 --fluid --fluid-pressure-to… | |
+ /bin/bash -c '\ | |
+ for f in diurnal.output*.txt; do \ | |
+ gnuplot -e "filename=\"$$f\"" $<_fluid.gp > $$f.png; \ | |
+ done' | |
+ | |
1d_fd_simple_shear: $(OBJ) $(HDR) | |
$(CC) $(LDFLAGS) $(OBJ) -o $@ | |
diff --git a/fluid.c b/fluid.c | |
t@@ -31,20 +31,29 @@ static double darcy_pressure_change_1d( | |
const double beta_f, | |
const double mu_f) | |
{ | |
- const double p = p_f_ghost_in[idx1g(i)]; | |
- double p_zn = p_f_ghost_in[idx1g(i-1)]; | |
- double p_zp = p_f_ghost_in[idx1g(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]; | |
- const double k_zn = k[i-1]; | |
- const double k_zp = k[i+1]; | |
- | |
- double div_k_grad_p = | |
- (2.*k_zp*k_/(k_zp + k_) * | |
- (p_zp - p)/dz | |
- - | |
- 2.*k_zn*k_/(k_zn + k_) * | |
- (p - p_zn)/dz)/dz; | |
+ 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); | |
+#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; | |
+#ifdef DEBUG | |
+ 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; | |
t@@ -104,7 +113,9 @@ int darcy_solver_1d( | |
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); | |
- /* puts(".. p_f_ghost after BC:"); print_array(sim->p_f_ghost, sim->nz… | |
+#ifdef DEBUG | |
+ 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) | |
t@@ -120,6 +131,11 @@ int darcy_solver_1d( | |
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]); | |
+#endif | |
+ | |
p_f = sim->p_f_ghost[idx1g(i)]; | |
p_f_ghost_out[idx1g(i)] = p_f | |
t@@ -133,20 +149,27 @@ int darcy_solver_1d( | |
} | |
r_norm_max = max(r_norm, sim->nz); | |
- /* puts(".. p_f_ghost_out:"); print_array(p_f_ghost_out, sim->nz+2); */ | |
+#ifdef DEBUG | |
+ 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); | |
- /* puts(".. p_f_ghost after update:"); print_array(sim->p_f_ghost, sim… | |
+#ifdef DEBUG | |
+ 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; /* Include top node in… | |
+ 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); | |
- /* printf(".. Solution converged after %d iterations\n", iter); */ | |
+#ifdef DEBUG | |
+ printf(".. Solution converged after %d iterations\n", iter); | |
+#endif | |
return 0; | |
} | |
} | |
diff --git a/main.c b/main.c | |
t@@ -42,7 +42,7 @@ static void usage(void) | |
" -t, --time VAL simulation start time [s]\n" | |
" -T, --time-end VAL simulation end time [s]\n" | |
" -D, --time-step VAL computational time step length [… | |
- " -I, --file-time VAL interval between output files [s… | |
+ " -I, --file-interval VAL interval between output files [s… | |
" -v, --version show version information\n" | |
" -h, --help show this message\n" | |
, __func__, PROGNAME); | |
t@@ -95,7 +95,7 @@ int main(int argc, char* argv[]) | |
{"time", required_argument, NULL, 't'}, | |
{"time-end", required_argument, NULL, 'T'}, | |
{"time-step", required_argument, NULL, 'D'}, | |
- {"file-time", required_argument, NULL, 'I'}, | |
+ {"file-interval", required_argument, NULL, 'I'}, | |
{NULL, 0, NULL, 0} | |
}; | |
t@@ -224,19 +224,24 @@ int main(int argc, char* argv[]) | |
if (sim.fluid) | |
hydrostatic_fluid_pressure_distribution(&sim); | |
- /* puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, sim… | |
- /* puts(""); */ | |
- /* puts(".. normal stress before iterations:"); print_array(sim.sigma_n, s… | |
- /* puts(""); */ | |
+#ifdef DEBUG | |
+ puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, sim.nz… | |
+ puts(""); | |
+ puts(".. normal stress before iterations:"); print_array(sim.sigma_n, sim.… | |
+ puts(""); | |
+#endif | |
+ write_output_file(&sim); | |
double filetimeclock = 0.0; | |
while (sim.t <= sim.t_end) { | |
if (sim.fluid) { | |
if (darcy_solver_1d(&sim, 10000, 1e-3)) | |
exit(1); | |
- /* puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2); */ | |
- /* puts(""); */ | |
+#ifdef DEBUG | |
+ puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2); | |
+ puts(""); | |
+#endif | |
} | |
compute_effective_stress(&sim); | |
diff --git a/simulation.c b/simulation.c | |
t@@ -397,7 +397,8 @@ void write_output_file(struct simulation* sim) | |
char outfile[200]; | |
FILE *fp; | |
- sprintf(outfile, "%s.output%05d.txt", sim->name, sim->n_file++); | |
+ sprintf(outfile, "%s.output%05d.t=%gs.txt", | |
+ sim->name, sim->n_file++, sim->t); | |
fp = fopen(outfile, "w"); | |
if (sim->fluid) |