commit edf97b34d7ba2dcd925192cd64479d16c4767b2a parent b2f8d047df70db4aa5593bfafe9bd594312e1481 Author: Anders Damsgaard <anders@adamsgaard.dk> Date: Fri, 28 Jun 2019 14:07:38 +0000 Merge branch 'debug_diffusion' into 'master' Debug diffusion See merge request admesg/1d_fd_simple_shear!2 Diffstat:
31 files changed, 516 insertions(+), 419 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml @@ -2,12 +2,10 @@ build-alpine: stage: build image: alpine:edge before_script: - - apk --no-cache add build-base gnuplot bash zsh valgrind + - apk --no-cache add build-base gnuplot valgrind script: - - make plots - - valgrind --error-exitcode=1 --leak-check=full ./1d_fd_simple_shear -h - - valgrind --error-exitcode=1 --leak-check=full ./1d_fd_simple_shear - - valgrind --error-exitcode=1 --leak-check=full ./1d_fd_simple_shear -F + - make -C examples + - make memtest artifacts: paths: - - 1d_fd_simple_shear.png + - examples/1d_fd_simple_shear.png diff --git a/1d_fd_simple_shear_henann_kamrin2016.h b/1d_fd_simple_shear_henann_kamrin2016.h @@ -1,36 +0,0 @@ -#ifndef ONED_FD_SIMPLE_SHEAR_ -#define ONED_FD_SIMPLE_SHEAR_ - -#include <math.h> -#include "arrays.h" -#include "simulation.h" - -#define DEG2RAD(x) (x*M_PI/180.0) - -/* Simulation settings */ -struct simulation init_sim(void) -{ - struct simulation sim; - - sim.G = 9.81; - - sim.P_wall = 10e3; /* larger normal stress deepens the deformational depth */ - sim.mu_wall = 0.40; - sim.v_x_bot = 0.0; - - sim.nz = 100; - - sim.A = 0.48; - sim.b = 0.9377; - sim.mu_s = 0.3819; - sim.phi = initval(0.38, sim.nz); - sim.d = 1e-3; - sim.rho_s = 2.485e3; - - sim.origo_z = 0.0; - sim.L_z = 20.0*sim.d; - - return sim; -} - -#endif diff --git a/Makefile b/Makefile @@ -3,16 +3,17 @@ LDFLAGS = -lm SRC = $(wildcard *.c) OBJ = $(patsubst %.c,%.o,$(SRC)) HDR = $(wildcard *.h) -BIN = 1d_fd_simple_shear +BIN = ./1d_fd_simple_shear PREFIX ?= /usr/local INSTALL ?= install STRIP ?= strip -.PHONY: default default: $(BIN) -.PHONY: install +$(BIN): $(OBJ) $(HDR) + $(CC) $(LDFLAGS) $(OBJ) -o $@ + install: $(BIN) $(STRIP) $(BIN) $(INSTALL) -m 0755 -d $(DESTDIR)$(PREFIX)/bin @@ -21,109 +22,17 @@ install: $(BIN) uninstall: $(RM) $(DESTDIR)$(PREFIX)/bin/$(BIN) -.PHONY: plots -plots: 1d_fd_simple_shear.png \ - 1d_fd_simple_shear_rheology.png \ - 1d_fd_simple_shear_rheology_kamb.png \ - 1d_fd_simple_shear_rheology_iverson.png \ - 1d_fd_simple_shear_rheology_tulaczyk.png \ - diurnal.timeseries.pdf - -diurnal.mp4: diurnal.output00000.txt.png - ffmpeg -i diurnal.output%05d.txt.png -y diurnal.mp4 - -diurnal.output00000.txt.png: 1d_fd_simple_shear_fluid.gp diurnal.output00000.txt - /bin/bash -c '\ - for f in diurnal.output*.txt; do \ - gnuplot -e "filename=\"$$f\"; p_min=\"0\"; p_max=\"100e3\"" $< > $$f.png; \ - done' - -diurnal.output00000.txt: 1d_fd_simple_shear - /usr/bin/env zsh -c '\ - ./$< --resolution 50 --length 2.0 --normal-stress 150e3 --fluid --fluid-permeability 2e-17 --fluid-pressure-top 50e3 --fluid-pressure-ampl 50e3 --fluid-pressure-freq $$(( 1.0/(3600*24) )) --time-step 1e-1 --file-interval $$(( 60*10 )) --time-end $$(( 3600*24*7 )) diurnal' - -diurnal.timeseries.txt: diurnal.output00000.txt - /bin/bash -c '\ - for f in diurnal.output*.txt; do \ - tail -n 1 "$$f" | cut -f2- >> $@; \ - done' - -diurnal.timeseries.pdf: diurnal.timeseries.gp diurnal.timeseries.txt - gnuplot $< > $@ - -diurnal.gif: diurnal.mp4 - convert diurnal.output*.txt.png \ - -delay 1 -loop 0 -fuzz 10% -layers Optimize $@ - -1d_fd_simple_shear: $(OBJ) $(HDR) - $(CC) $(LDFLAGS) $(OBJ) -o $@ - -1d_fd_simple_shear.png: 1d_fd_simple_shear 1d_fd_simple_shear.gp - /bin/bash -c '\ - for P in 10 20 40 60 80 120; do \ - ./$< -o 0.03 -L 0.64 -P $${P}e3 -N > $<_P$${P}kPa.txt; \ - done' - gnuplot $<.gp > $@ - -1d_fd_simple_shear_rheology.png: 1d_fd_simple_shear 1d_fd_simple_shear_rheology.gp - /bin/bash -c '\ - for b in $$(printf "0.01\n0.10\n"; seq 0.20 0.20 1.00); do \ - out="$<_rheology_b$${b}.txt"; \ - rm -f "$$out"; \ - for t in $$(seq 0.01 0.01 0.8); do \ - printf "$$t\t" >> "$$out"; \ - ./$< -P 20e3 --stress-ratio $$t -b $$b | \ - tail -n 1 | cut -f2 >> "$$out"; \ - done; done' - gnuplot $<_rheology.gp > $@ - -# 1 bar is equal to 100 kPa -1d_fd_simple_shear_rheology_kamb.png: 1d_fd_simple_shear 1d_fd_simple_shear_rheology_kamb.gp - /bin/bash -c '\ - for b in 0.01 0.10 0.20 0.40 0.94; do \ - out="$<_rheology_b$${b}_kamb.txt"; \ - rm -f "$$out"; \ - for t in $$(seq 0.01 0.01 2.0); do \ - printf "$$t\t" >> "$$out"; \ - ./$< -f 1.1 -P 1.7e3 --stress-ratio $$t -b $$b | \ - tail -n 1 | cut -f2 >> "$$out"; \ - done; done' - gnuplot $<_rheology_kamb.gp > $@ - -# shear-strain rate from 10^1 to 10^6 m/a -# friction around 0.55 -1d_fd_simple_shear_rheology_iverson.png: 1d_fd_simple_shear 1d_fd_simple_shear_rheology_iverson.gp - /bin/bash -c '\ - for b in $$(printf "0.01\n0.10\n"; seq 0.20 0.20 0.90) 0.94; do \ - out="$<_rheology_b$${b}_iverson.txt"; \ - rm -f "$$out"; \ - for t in $$(seq 0.0001 0.002 1.0); do \ - printf "$$t\t" >> "$$out"; \ - ./$< -f 0.55 -P 100e3 -L 1.0 --stress-ratio $$t -b $$b | \ - tail -n 1 | cut -f2 >> "$$out"; \ - done; done' - gnuplot $<_rheology_iverson.gp > $@ - -# shear velocity rate from 0.1 m/h to 100 m/h -1d_fd_simple_shear_rheology_tulaczyk.png: 1d_fd_simple_shear 1d_fd_simple_shear_rheology_tulaczyk.gp - /bin/bash -c '\ - for b in $$(printf "0.01\n0.10\n"; seq 0.20 0.20 1.00) 0.94; do \ - out="$<_rheology_b$${b}_tulaczyk.txt"; \ - rm -f "$$out"; \ - for t in $$(seq 0.1 0.002 0.9); do \ - printf "$$t\t" >> "$$out"; \ - ./$< -f 0.5 -P 10e3 --stress-ratio $$t -b $$b | \ - tail -n 1 | cut -f2 >> "$$out"; \ - done; done' - gnuplot $<_rheology_tulaczyk.gp > $@ - -.PHONY: watch watch: echo $(SRC) $(HDR) | tr ' ' '\n' | entr -s 'make && ./1d_fd_simple_shear' +memtest: $(BIN) + valgrind --error-exitcode=1 --leak-check=full $(BIN) -h + valgrind --error-exitcode=1 --leak-check=full $(BIN) + valgrind --error-exitcode=1 --leak-check=full $(BIN) -F + clean: $(RM) *.txt $(RM) *.o $(RM) 1d_fd_simple_shear - $(RM) 1d_fd_simple_shear.png - $(RM) 1d_fd_simple_shear_rheology*.png + +.PHONY: default install uninstall watch memtest clean diff --git a/README.md b/README.md @@ -59,7 +59,7 @@ probably be improved further. | Discrete-element model | Continuum Model | | ----------------------- | --------------- | | Damsgaard et al. 2013 | | -|  |  | +|  |  | ### Stress and strain rate The rheology is of Bingham type, where no deformation occurs beneath the @@ -71,15 +71,15 @@ have *b* = 0.94 ([Forterre and Pouliquen | Real material (laboratory or field study) | Continuum Model | | -------------------------------------------- | --------------- | | Upstream-B ([Kamb 1991](https://doi.org/10.1029/91jb00946)): | | -|  |  | +|  |  | | Various subglacial tills ([Iverson 2010](https://doi.org/10.3189/002214311796406220)): | | -|  |  | +|  |  | | Whillans Ice Plain ([Tulaczyk 2006](https://doi.org/10.3189/172756506781828601)): | | -|  |  | +|  |  | ### Variable water pressure The model is expanded from the Henann and Kamrin 2013 model by including a diffusive porewater pressure parameterization. Below is an example of diurnal water-pressure variations that gradually propagate into the bed. - + diff --git a/arrays.c b/arrays.c @@ -6,158 +6,214 @@ /* Translate a i,j,k index in grid with dimensions nx, ny, nz into a linear * index */ -unsigned int idx3(const unsigned int i, - const unsigned int j, - const unsigned int k, - const unsigned int nx, - const unsigned int ny) +unsigned int +idx3(const unsigned int i, + const unsigned int j, + const unsigned int k, + const unsigned int nx, + const unsigned int ny) { return i + nx*j + nx*ny*k; } /* Translate a i,j,k index in grid with dimensions nx, ny, nz and a padding of * single ghost nodes into a linear index */ -unsigned int idx3g(const unsigned int i, - const unsigned int j, - const unsigned int k, - const unsigned int nx, - const unsigned int ny) +unsigned int +idx3g(const unsigned int i, + const unsigned int j, + const unsigned int k, + const unsigned int nx, + const unsigned int ny) { return i+1 + (nx+2)*(j+1) + (nx+2)*(ny+2)*(k+1); } /* Translate a i,j index in grid with dimensions nx, ny into a linear index */ -unsigned int idx2(const unsigned int i, - const unsigned int j, - const unsigned int nx) +unsigned int +idx2(const unsigned int i, const unsigned int j, const unsigned int nx) { return i + nx*j; } /* Translate a i,j index in grid with dimensions nx, ny and a padding of single * ghost nodes into a linear index */ -unsigned int idx2g(const unsigned int i, - const unsigned int j, - const unsigned int nx) +unsigned int +idx2g(const unsigned int i, const unsigned int j, const unsigned int nx) { return i+1 + (nx+2)*(j+1); } /* Translate a i index in grid with a padding of single into a linear index */ -unsigned int idx1g(const unsigned int i) +unsigned int +idx1g(const unsigned int i) { return i+1; } /* Return an array of `n` linearly spaced values in the range [lower; upper] */ -double* linspace(const double lower, const double upper, const int n) +double* +linspace(const double lower, const double upper, const int n) { - double *x = malloc(n*sizeof(double)); - double dx = (upper - lower)/(double)(n-1); - for (int i=0; i<n; ++i) + int i; + double *x; + double dx; + + x = malloc(n*sizeof(double)); + dx = (upper - lower)/(double)(n-1); + for (i=0; i<n; ++i) x[i] = lower + dx*i; return x; } +/* Return an array of `n-1` values with the intervals between `x` values */ +double* +spacing(const double* x, const int n) +{ + int i; + double *dx; + + dx = malloc((n-1)*sizeof(double)); + for (i=0; i<n-1; ++i) + dx[i] = x[i+1] - x[i]; + return dx; +} + /* Return an array of `n` values with the value 0.0 */ -double* zeros(const int n) +double* +zeros(const int n) { - double *x = malloc(n*sizeof(double)); - for (int i=0; i<n; ++i) + int i; + double *x; + + x = malloc(n*sizeof(double)); + for (i=0; i<n; ++i) x[i] = 0.0; return x; } /* Return an array of `n` values with the value 1.0 */ -double* ones(const int n) +double* +ones(const int n) { - double *x = malloc(n*sizeof(double)); - for (int i=0; i<n; ++i) + int i; + double *x; + + x = malloc(n*sizeof(double)); + for (i=0; i<n; ++i) x[i] = 1.0; return x; } /* Return an array of `n` values with a specified value */ -double* initval(const double value, const int n) +double* +initval(const double value, const int n) { - double *x = malloc(n*sizeof(double)); - for (int i=0; i<n; ++i) + int i; + double *x; + + x = malloc(n*sizeof(double)); + for (i=0; i<n; ++i) x[i] = value; return x; } /* Return an array of `n` uninitialized values */ -double* empty(const int n) +double* +empty(const int n) { return malloc(n*sizeof(double)); } /* Return largest value in array `a` with size `n` */ -double max(const double* a, const int n) +double +max(const double* a, const int n) { - double maxval = -INFINITY; - for (int i=0; i<n; ++i) + int i; + double maxval; + + maxval = -INFINITY; + for (i=0; i<n; ++i) if (a[i] > maxval) maxval = a[i]; return maxval; } /* Return smallest value in array `a` with size `n` */ -double min(const double* a, const int n) +double +min(const double* a, const int n) { - double minval = +INFINITY; - for (int i=0; i<n; ++i) + int i; + double minval; + + minval = +INFINITY; + for (i=0; i<n; ++i) if (a[i] < minval) minval = a[i]; return minval; } -void print_array(const double* a, const int n) +void +print_array(const double* a, const int n) { - for (int i=0; i<n; ++i) + int i; + for (i=0; i<n; ++i) printf("%.17g\n", a[i]); } -void print_arrays(const double* a, const double* b, const int n) +void +print_arrays(const double* a, const double* b, const int n) { - for (int i=0; i<n; ++i) + int i; + for (i=0; i<n; ++i) printf("%.17g\t%.17g\n", a[i], b[i]); } -void print_arrays_2nd_normalized(const double* a, const double* b, const int n) +void +print_arrays_2nd_normalized(const double* a, const double* b, const int n) { - double max_b = max(b, n); - for (int i=0; i<n; ++i) + int i; + double max_b; + + max_b = max(b, n); + for (i=0; i<n; ++i) printf("%.17g\t%.17g\n", a[i], b[i]/max_b); } -void print_three_arrays(const double* a, - const double* b, - const double* c, - const int n) +void +print_three_arrays(const double* a, + const double* b, + const double* c, + const int n) { - for (int i=0; i<n; ++i) + int i; + for (i=0; i<n; ++i) printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); } -void fprint_arrays(FILE* fp, const double* a, const double* b, const int n) +void +fprint_arrays(FILE* fp, const double* a, const double* b, const int n) { - for (int i=0; i<n; ++i) + int i; + for (i=0; i<n; ++i) fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]); } -void fprint_three_arrays(FILE* fp, - const double* a, - const double* b, - const double* c, - const int n) +void +fprint_three_arrays(FILE* fp, + const double* a, + const double* b, + const double* c, + const int n) { - for (int i=0; i<n; ++i) + int i; + for (i=0; i<n; ++i) fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); } -void copy_values(const double* in, double* out, const int n) +void +copy_values(const double* in, double* out, const int n) { - for (int i=0; i<n; ++i) + int i; + for (i=0; i<n; ++i) out[i] = in[i]; } diff --git a/arrays.h b/arrays.h @@ -19,6 +19,7 @@ unsigned int idx2g( unsigned int idx1g(const unsigned int i); +double* spacing(const double* x, const int n); double* linspace(const double lower, const double upper, const int n); double* zeros(const int n); double* ones(const int n); diff --git a/diurnal.gif b/diurnal.gif Binary files differ. diff --git a/damsgaard2013-fig8.png b/doc/damsgaard2013-fig8.png Binary files differ. diff --git a/iverson2010-fig2a.png b/doc/iverson2010-fig2a.png Binary files differ. diff --git a/iverson2010-fig2b.png b/doc/iverson2010-fig2b.png Binary files differ. diff --git a/kamb1991-fig1.png b/doc/kamb1991-fig1.png Binary files differ. diff --git a/tulaczyk2006-fig1.png b/doc/tulaczyk2006-fig1.png Binary files differ. diff --git a/1d_fd_simple_shear.gp b/examples/1d_fd_simple_shear.gp diff --git a/1d_fd_simple_shear.png b/examples/1d_fd_simple_shear.png Binary files differ. diff --git a/1d_fd_simple_shear_fluid.gp b/examples/1d_fd_simple_shear_fluid.gp diff --git a/1d_fd_simple_shear_rheology.gp b/examples/1d_fd_simple_shear_rheology.gp diff --git a/1d_fd_simple_shear_rheology.png b/examples/1d_fd_simple_shear_rheology.png Binary files differ. diff --git a/1d_fd_simple_shear_rheology_iverson.gp b/examples/1d_fd_simple_shear_rheology_iverson.gp diff --git a/1d_fd_simple_shear_rheology_iverson.png b/examples/1d_fd_simple_shear_rheology_iverson.png Binary files differ. diff --git a/1d_fd_simple_shear_rheology_kamb.gp b/examples/1d_fd_simple_shear_rheology_kamb.gp diff --git a/1d_fd_simple_shear_rheology_kamb.png b/examples/1d_fd_simple_shear_rheology_kamb.png Binary files differ. diff --git a/1d_fd_simple_shear_rheology_tulaczyk.gp b/examples/1d_fd_simple_shear_rheology_tulaczyk.gp diff --git a/1d_fd_simple_shear_rheology_tulaczyk.png b/examples/1d_fd_simple_shear_rheology_tulaczyk.png Binary files differ. diff --git a/BlueSeq.plt b/examples/BlueSeq.plt diff --git a/examples/Makefile b/examples/Makefile @@ -0,0 +1,105 @@ +BIN = ../1d_fd_simple_shear + +.PHONY: plots +plots: 1d_fd_simple_shear.png \ + 1d_fd_simple_shear_rheology.png \ + 1d_fd_simple_shear_rheology_kamb.png \ + 1d_fd_simple_shear_rheology_iverson.png \ + 1d_fd_simple_shear_rheology_tulaczyk.png \ + diurnal.timeseries.pdf + +diurnal.mp4: diurnal.output00000.txt.png + ffmpeg -i diurnal.output%05d.txt.png -y diurnal.mp4 + +diurnal.output00000.txt.png: 1d_fd_simple_shear_fluid.gp diurnal.output00000.txt + /bin/sh -c '\ + for f in diurnal.output*.txt; do \ + gnuplot -e "filename=\"$$f\"; p_min=\"0\"; p_max=\"100e3\"" $< > $$f.png; \ + done' + +diurnal.output00000.txt: $(BIN) + /bin/sh -c '\ + ./$< --resolution 50 --length 2.0 --normal-stress 150e3 --fluid --fluid-permeability 2e-17 --fluid-pressure-top 50e3 --fluid-pressure-ampl 50e3 --fluid-pressure-freq $$( echo "1.0/(3600*24)" | bc -l ) --file-interval $$( echo "60*10" | bc -l ) --time-end $$( echo "3600*24*7" | bc -l ) diurnal' + +diurnal.timeseries.txt: diurnal.output00000.txt + /bin/sh -c '\ + for f in diurnal.output*.txt; do \ + tail -n 1 "$$f" | cut -f2- >> $@; \ + done' + +diurnal.timeseries.pdf: diurnal.timeseries.gp diurnal.timeseries.txt + gnuplot $< > $@ + +diurnal.gif: diurnal.mp4 + convert diurnal.output*.txt.png \ + -delay 1 -loop 0 -fuzz 10% -layers Optimize $@ + +1d_fd_simple_shear.png: $(BIN) 1d_fd_simple_shear.gp + /bin/sh -c '\ + for P in 10 20 40 60 80 120; do \ + ./$< -o 0.03 -L 0.64 -P $${P}e3 -N > 1d_fd_simple_shear_P$${P}kPa.txt; \ + done' + gnuplot 1d_fd_simple_shear.gp > $@ + +1d_fd_simple_shear_rheology.png: $(BIN) 1d_fd_simple_shear_rheology.gp + /bin/sh -c '\ + for b in $$(printf "0.01\n0.10\n"; seq 0.20 0.20 1.00); do \ + out="1d_fd_simple_shear_rheology_b$${b}.txt"; \ + rm -f "$$out"; \ + for t in $$(seq 0.01 0.01 0.8); do \ + printf "$$t\t" >> "$$out"; \ + ./$< -P 20e3 --stress-ratio $$t -b $$b | \ + tail -n 1 | cut -f2 >> "$$out"; \ + done; done' + gnuplot 1d_fd_simple_shear_rheology.gp > $@ + +# 1 bar is equal to 100 kPa +1d_fd_simple_shear_rheology_kamb.png: $(BIN) 1d_fd_simple_shear_rheology_kamb.gp + /bin/sh -c '\ + for b in 0.01 0.10 0.20 0.40 0.94; do \ + out="1d_fd_simple_shear_rheology_b$${b}_kamb.txt"; \ + rm -f "$$out"; \ + for t in $$(seq 0.01 0.01 2.0); do \ + printf "$$t\t" >> "$$out"; \ + ./$< -f 1.1 -P 1.7e3 --stress-ratio $$t -b $$b | \ + tail -n 1 | cut -f2 >> "$$out"; \ + done; done' + gnuplot 1d_fd_simple_shear_rheology_kamb.gp > $@ + +# shear-strain rate from 10^1 to 10^6 m/a +# friction around 0.55 +1d_fd_simple_shear_rheology_iverson.png: $(BIN) 1d_fd_simple_shear_rheology_iverson.gp + /bin/sh -c '\ + for b in $$(printf "0.01\n0.10\n"; seq 0.20 0.20 0.90) 0.94; do \ + out="1d_fd_simple_shear_rheology_b$${b}_iverson.txt"; \ + rm -f "$$out"; \ + for t in $$(seq 0.0001 0.002 1.0); do \ + printf "$$t\t" >> "$$out"; \ + ./$< -f 0.55 -P 100e3 -L 1.0 --stress-ratio $$t -b $$b | \ + tail -n 1 | cut -f2 >> "$$out"; \ + done; done' + gnuplot 1d_fd_simple_shear_rheology_iverson.gp > $@ + +# shear velocity rate from 0.1 m/h to 100 m/h +1d_fd_simple_shear_rheology_tulaczyk.png: $(BIN) 1d_fd_simple_shear_rheology_tulaczyk.gp + /bin/sh -c '\ + for b in $$(printf "0.01\n0.10\n"; seq 0.20 0.20 1.00) 0.94; do \ + out="1d_fd_simple_shear_rheology_b$${b}_tulaczyk.txt"; \ + rm -f "$$out"; \ + for t in $$(seq 0.1 0.002 0.9); do \ + printf "$$t\t" >> "$$out"; \ + ./$< -f 0.5 -P 10e3 --stress-ratio $$t -b $$b | \ + tail -n 1 | cut -f2 >> "$$out"; \ + done; done' + gnuplot 1d_fd_simple_shear_rheology_tulaczyk.gp > $@ + +$(BIN): + make -C .. + +clean: + $(RM) *.txt + $(RM) *.png + $(RM) *.pdf + $(RM) *.o + $(RM) 1d_fd_simple_shear.png + $(RM) 1d_fd_simple_shear_rheology*.png diff --git a/examples/diurnal.gif b/examples/diurnal.gif Binary files differ. diff --git a/diurnal.timeseries.gp b/examples/diurnal.timeseries.gp diff --git a/fluid.c b/fluid.c @@ -3,39 +3,81 @@ #include "simulation.h" #include "arrays.h" -void hydrostatic_fluid_pressure_distribution(struct simulation* sim) +/* #define DEBUG true */ + +void +hydrostatic_fluid_pressure_distribution(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) + int i; + for (i=0; i<sim->nz; ++i) sim->p_f_ghost[idx1g(i)] = sim->p_f_top + sim->phi[i]*sim->rho_f*sim->G* (sim->L_z - sim->z[i]); } -static double sine_wave(const double time, - const double amplitude, - const double frequency, - const double phase, - const double base_value) +/* Determines the largest time step for the current simulation state. Note + * that the time step should be recalculated if cell sizes or diffusivities + * (i.e., permeabilities, porosities, viscosities, or compressibilities) + * change. The safety factor should be in ]0;1] */ +void +set_largest_fluid_timestep(struct simulation* sim, const double safety) +{ + int i; + double dx_min, diff, diff_max; + double *dx; + + dx = spacing(sim->z, sim->nz); + dx_min = INFINITY; + for (i=0; i<sim->nz-1; ++i) { + if (dx[i] < 0.0) { + fprintf(stderr, "error: cell spacing negative (%g) in cell %d\n", + dx[i], i); + free(dx); + exit(10); + } + if (dx[i] < dx_min) dx_min = dx[i]; + } + free(dx); + + diff_max = -INFINITY; + for (i=0; i<sim->nz; ++i) { + diff = sim->k[i]/(sim->phi[i]*sim->beta_f*sim->mu_f); + if (diff > diff_max) diff_max = diff; + } + + sim->dt = safety*0.5*dx_min*dx_min/diff_max; + if (sim->file_dt*0.5 < sim->dt) + sim->dt = sim->file_dt; +} + +static double +sine_wave(const double time, + const double amplitude, + const double frequency, + const double phase, + const double base_value) { return amplitude*sin(2.0*PI*frequency*time + phase) + base_value; } -static double darcy_pressure_change_1d(const int i, - const int nz, - const double* p_f_ghost_in, - const double* phi, - const double* k, - const double dz, - const double dt, - const double beta_f, - const double mu_f) +static double +darcy_pressure_change_1d(const int i, + const int nz, + const double* p_f_ghost_in, + const double* phi, + const double* k, + const double dz, + const double dt, + const double beta_f, + const double mu_f) { - 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)]; + double p, p_zn, p_zp, k_, div_k_grad_p, k_zn, k_zp; - const double k_ = k[i]; - double k_zn, k_zp; + p = p_f_ghost_in[idx1g(i)]; + p_zn = p_f_ghost_in[idx1g(i-1)]; + p_zp = p_f_ghost_in[idx1g(i+1)]; + + k_ = k[i]; 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 @@ -45,67 +87,71 @@ static double darcy_pressure_change_1d(const int i, 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; + 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; + return dt/(beta_f*phi[i]*mu_f)*div_k_grad_p; } -int darcy_solver_1d(struct simulation* sim, - const int max_iter, - const double rel_tol) +int +darcy_solver_1d(struct simulation* sim, + const int max_iter, + const double rel_tol) { - - /* compute explicit solution to pressure change */ - double* dp_f_expl = zeros(sim->nz); - for (int i=0; i<sim->nz; ++i) - dp_f_expl[i] = darcy_pressure_change_1d(i, - sim->nz, - sim->p_f_ghost, - sim->phi, - sim->k, - sim->dz, - sim->dt, - sim->beta_f, - sim->mu_f); + int i, iter; + double epsilon, theta, p_f, p_f_top, r_norm_max; + double *dp_f_expl, *dp_f_impl, *p_f_ghost_out, *r_norm; /* choose integration method, parameter in [0.0; 1.0] * epsilon = 0.0: explicit * epsilon = 0.5: Crank-Nicolson * epsilon = 1.0: implicit */ - const double epsilon = 0.5; + epsilon = 0.5; /* choose relaxation factor, parameter in ]0.0; 1.0] * theta in ]0.0; 1.0]: underrelaxation * theta = 1.0: Gauss-Seidel * theta > 1.0: overrelaxation */ - const double theta = 0.05; - /* const double theta = 1.7; */ + /* TODO: values other than 1.0 do not work! */ + theta = 1.0; - double p_f; + p_f_top = sine_wave(sim->t, + sim->p_f_mod_ampl, + sim->p_f_mod_freq, + sim->p_f_mod_phase, + sim->p_f_top); - /* compute implicit solution to pressure change */ - int iter; - double* dp_f_impl = zeros(sim->nz); - double* p_f_ghost_out = zeros(sim->nz+2); - double* r_norm = zeros(sim->nz); - double r_norm_max = NAN; - double p_f_top = sine_wave( - sim->t, - sim->p_f_mod_ampl, - sim->p_f_mod_freq, - sim->p_f_mod_phase, - sim->p_f_top); + /* set fluid BCs (1 of 2) */ + 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); + + /* compute explicit solution to pressure change */ + dp_f_expl = zeros(sim->nz); + for (i=0; i<sim->nz; ++i) + dp_f_expl[i] = darcy_pressure_change_1d(i, + sim->nz, + sim->p_f_ghost, + sim->phi, + sim->k, + sim->dz, + sim->dt, + sim->beta_f, + sim->mu_f); + /* compute implicit solution to pressure change */ + dp_f_impl = zeros(sim->nz); + p_f_ghost_out = zeros(sim->nz+2); + r_norm = zeros(sim->nz); for (iter=0; iter<max_iter; ++iter) { + /* set fluid BCs (2 of 2) */ 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); @@ -114,7 +160,7 @@ int darcy_solver_1d(struct simulation* sim, #endif /* for (int i=0; i<sim->nz; ++i) */ - for (int i=0; i<sim->nz-1; ++i) + for (i=0; i<sim->nz-1; ++i) dp_f_impl[i] = darcy_pressure_change_1d(i, sim->nz, sim->p_f_ghost, @@ -124,8 +170,7 @@ int darcy_solver_1d(struct simulation* sim, sim->dt, sim->beta_f, sim->mu_f); - /* for (int i=0; i<sim->nz; ++i) { */ - for (int i=0; i<sim->nz-1; ++i) { + for (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]); @@ -155,7 +200,7 @@ int darcy_solver_1d(struct simulation* sim, print_array(sim->p_f_ghost, sim->nz+2); #endif - if (r_norm_max <= rel_tol) { + 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; /* top node in BC */ set_bc_neumann(sim->p_f_ghost, sim->nz, -1); @@ -167,7 +212,7 @@ int darcy_solver_1d(struct simulation* sim, printf(".. Solution converged after %d iterations\n", iter); #endif return 0; - } + } } free(dp_f_expl); diff --git a/fluid.h b/fluid.h @@ -5,6 +5,8 @@ void hydrostatic_fluid_pressure_distribution(struct simulation* sim); +void set_largest_fluid_timestep(struct simulation* sim, const double safety); + int darcy_solver_1d( struct simulation* sim, const int max_iter, diff --git a/main.c b/main.c @@ -6,18 +6,21 @@ #include "simulation.h" #include "fluid.h" -#define VERSION "0.2" +#define VERSION "0.3.0" #define PROGNAME "1d_fd_simple_shear" #include "parameter_defaults.h" -static void usage(void) +static void +usage(void) { struct simulation sim = init_sim(); printf("%s: %s [OPTIONS] [NAME]\n" "runs a simulation and outputs state in files prefixed with NAME.\n" "If NAME is not specified, the default value '%s' is used.\n" - "optional arguments:\n" + "\nOptional arguments:\n" + " -v, --version show version information\n" + " -h, --help show this message\n" " -N, --normalize normalize output velocity\n" " -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n" " -P, --normal-stress VAL normal stress on top [Pa] (default %g)\n" @@ -32,6 +35,7 @@ static void usage(void) " -n, --resolution VAL number of cells in domain [-] (default %d)\n" " -o, --origo VAL coordinate system origo [m] (default %g)\n" " -L, --length VAL domain length [m] (default %g)\n" + "\nOptional arguments only relevant with transient (fluid) simulation:\n" " -F, --fluid enable pore fluid computations\n" " -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default %g)\n" " -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default %g)\n" @@ -43,10 +47,7 @@ static void usage(void) " -H, --fluid-pressure-phase VAL fluid pressure at +z edge [Pa] (default %g)\n" " -t, --time VAL simulation start time [s] (default %g)\n" " -T, --time-end VAL simulation end time [s] (default %g)\n" - " -D, --time-step VAL computational time step length [s] (default %g)\n" - " -I, --file-interval VAL interval between output files [s] (default %g)\n" - " -v, --version show version information\n" - " -h, --help show this message\n", + " -I, --file-interval VAL interval between output files [s] (default %g)\n", __func__, PROGNAME, sim.name, sim.G, @@ -72,13 +73,13 @@ static void usage(void) sim.p_f_mod_phase, sim.t, sim.t_end, - sim.dt, sim.file_dt); free(sim.phi); free(sim.k); } -static void version(void) +static void +version(void) { printf("%s v%s\n" "Licensed under the GNU Public License, v3+\n" @@ -87,16 +88,21 @@ static void version(void) , PROGNAME, VERSION); } -int main(int argc, char* argv[]) +int +main(int argc, char* argv[]) { + int normalize, opt, i; + struct simulation sim; + const char* optstring; + unsigned long iter; + double new_phi, new_k, filetimeclock, max_v_x; /* load with default values */ - struct simulation sim = init_sim(); + sim = init_sim(); - int normalize = 0; + normalize = 0; - int opt; - const char* optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T:D:I:"; + optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T:D:I:"; const struct option longopts[] = { {"help", no_argument, NULL, 'h'}, {"version", no_argument, NULL, 'v'}, @@ -125,13 +131,12 @@ int main(int argc, char* argv[]) {"fluid-pressure-phase", required_argument, NULL, 'H'}, {"time", required_argument, NULL, 't'}, {"time-end", required_argument, NULL, 'T'}, - {"time-step", required_argument, NULL, 'D'}, {"file-interval", required_argument, NULL, 'I'}, {NULL, 0, NULL, 0} }; - double new_phi = sim.phi[0]; - double new_k = sim.k[0]; + new_phi = sim.phi[0]; + new_k = sim.k[0]; while ((opt = getopt_long(argc, argv, optstring, longopts, NULL)) != -1) { switch (opt) { case -1: /* no more arguments */ @@ -223,13 +228,9 @@ int main(int argc, char* argv[]) case 'T': sim.t_end = atof(optarg); break; - case 'D': - sim.dt = atof(optarg); - break; case 'I': sim.file_dt = atof(optarg); break; - default: fprintf(stderr, "%s: invalid option -- %c\n", argv[0], opt); fprintf(stderr, "Try `%s --help` for more information\n", @@ -237,7 +238,7 @@ int main(int argc, char* argv[]) return -2; } } - for (int i=optind; i<argc; ++i) { + for (i=optind; i<argc; ++i) { if (i>optind) { fprintf(stderr, "error: more than one simulation name specified\n"); @@ -249,50 +250,34 @@ int main(int argc, char* argv[]) prepare_arrays(&sim); if (!isnan(new_phi)) - for (int i=0; i<sim.nz; ++i) + for (i=0; i<sim.nz; ++i) sim.phi[i] = new_phi; if (!isnan(new_k)) - for (int i=0; i<sim.nz; ++i) + for (i=0; i<sim.nz; ++i) sim.k[i] = new_k; lithostatic_pressure_distribution(&sim); - if (sim.fluid) + if (sim.fluid) { hydrostatic_fluid_pressure_distribution(&sim); + set_largest_fluid_timestep(&sim, 0.5); + } -#ifdef DEBUG - puts(".. p_f_ghost before iterations:"); print_array(sim.p_f_ghost, sim.nz+2); - puts(""); - puts(".. normal stress before iterations:"); print_array(sim.sigma_n, sim.nz); - puts(""); -#endif + check_simulation_parameters(&sim); - double filetimeclock = 0.0; - unsigned long iter = 0; + filetimeclock = 0.0; + iter = 0; while (sim.t <= sim.t_end) { if (sim.fluid) { if (darcy_solver_1d(&sim, 10000, 1e-5)) exit(1); -#ifdef DEBUG - puts(".. p_f_ghost:"); print_array(sim.p_f_ghost, sim.nz+2); - puts(""); -#endif } compute_effective_stress(&sim); compute_friction(&sim); compute_cooperativity_length(&sim); - if (iter == 0) - check_simulation_parameters(&sim); - -#ifdef DEBUG - puts("\n## Before solver"); - puts(".. sigma_n_eff:"); print_array(sim.sigma_n_eff, sim.nz); - puts(".. mu:"); print_array(sim.mu, sim.nz); -#endif - if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5)) exit(1); @@ -310,13 +295,13 @@ int main(int argc, char* argv[]) } if (normalize) { - double max_v_x = max(sim.v_x, sim.nz); - for (int i=0; i<sim.nz; ++i) + max_v_x = max(sim.v_x, sim.nz); + for (i=0; i<sim.nz; ++i) sim.v_x[i] /= max_v_x; } if (sim.fluid) - for (int i=0; i<sim.nz; ++i) + for (i=0; i<sim.nz; ++i) printf("%.17g\t%.17g\t%.17g\t%.17g\n", sim.z[i], sim.v_x[i], diff --git a/simulation.c b/simulation.c @@ -4,7 +4,8 @@ #include "arrays.h" #include "simulation.h" -void prepare_arrays(struct simulation* sim) +void +prepare_arrays(struct simulation* sim) { sim->z = linspace(sim->origo_z, /* spatial coordinates */ sim->origo_z + sim->L_z, @@ -24,7 +25,8 @@ void prepare_arrays(struct simulation* sim) sim->g_ghost = zeros(sim->nz+2); /* fluidity with ghost nodes */ } -void free_arrays(struct simulation* sim) +void +free_arrays(struct simulation* sim) { free(sim->z); free(sim->mu); @@ -39,19 +41,18 @@ void free_arrays(struct simulation* sim) free(sim->g_ghost); } -static void warn_parameter_value( - const char message[], - const double value, - int* return_status) +static void +warn_parameter_value(const char message[], + const double value, + int* return_status) { fprintf(stderr, "check_simulation_parameters: %s (%.17g)\n", message, value); *return_status = 1; } -static void check_float(const char name[], - const double value, - int* return_status) +static void +check_float(const char name[], const double value, int* return_status) { if (isnan(value)) { warn_parameter_value("%s is NaN", value, return_status); @@ -62,9 +63,12 @@ static void check_float(const char name[], } } -void check_simulation_parameters(const struct simulation* sim) +void +check_simulation_parameters(const struct simulation* sim) { - int return_status = 0; + int return_status; + + return_status = 0; check_float("sim.G", sim->G, &return_status); if (sim->G < 0.0) @@ -160,31 +164,30 @@ void check_simulation_parameters(const struct simulation* sim) warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive", sim->p_f_mod_ampl, &return_status); - check_float("sim.p_f_mod_freq", sim->p_f_mod_freq, &return_status); - if (sim->p_f_mod_freq < 0.0) - warn_parameter_value("sim.p_f_mod_freq is not a zero or positive", - sim->p_f_mod_freq, &return_status); - - check_float("sim.beta_f", sim->beta_f, &return_status); - if (sim->beta_f <= 0.0) - warn_parameter_value("sim.beta_f is not positive", - sim->beta_f, &return_status); - - check_float("sim.mu_f", sim->mu_f, &return_status); - if (sim->mu_f <= 0.0) - warn_parameter_value("sim.mu_f is not positive", - sim->mu_f, &return_status); - - check_float("sim.rho_f", sim->rho_f, &return_status); - if (sim->rho_f <= 0.0) - warn_parameter_value("sim.rho_f is not positive", - sim->rho_f, &return_status); - - check_float("sim.k[0]", sim->k[0], &return_status); - if (sim->k[0] <= 0.0) - warn_parameter_value("sim.k[0] is not positive", - sim->k[0], &return_status); - + check_float("sim.p_f_mod_freq", sim->p_f_mod_freq, &return_status); + if (sim->p_f_mod_freq < 0.0) + warn_parameter_value("sim.p_f_mod_freq is not a zero or positive", + sim->p_f_mod_freq, &return_status); + + check_float("sim.beta_f", sim->beta_f, &return_status); + if (sim->beta_f <= 0.0) + warn_parameter_value("sim.beta_f is not positive", + sim->beta_f, &return_status); + + check_float("sim.mu_f", sim->mu_f, &return_status); + if (sim->mu_f <= 0.0) + warn_parameter_value("sim.mu_f is not positive", + sim->mu_f, &return_status); + + check_float("sim.rho_f", sim->rho_f, &return_status); + if (sim->rho_f <= 0.0) + warn_parameter_value("sim.rho_f is not positive", + sim->rho_f, &return_status); + + check_float("sim.k[0]", sim->k[0], &return_status); + if (sim->k[0] <= 0.0) + warn_parameter_value("sim.k[0] is not positive", + sim->k[0], &return_status); } if (return_status != 0) { @@ -193,83 +196,96 @@ void check_simulation_parameters(const struct simulation* sim) } } -void lithostatic_pressure_distribution(struct simulation* sim) +void +lithostatic_pressure_distribution(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) + int i; + for (i=0; i<sim->nz; ++i) sim->sigma_n[i] = sim->P_wall + (1.0 - sim->phi[i])*sim->rho_s*sim->G* (sim->L_z - sim->z[i]); } -void compute_friction(struct simulation* sim) +void +compute_friction(struct simulation* sim) { + int i; if (sim->fluid) - for (int i=0; i<sim->nz; ++i) + for (i=0; i<sim->nz; ++i) sim->mu[i] = sim->mu_wall/ (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top)); - else - for (int i=0; i<sim->nz; ++i) + for (i=0; i<sim->nz; ++i) sim->mu[i] = sim->mu_wall / (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G* (sim->L_z - sim->z[i])/sim->P_wall); } -double shear_strain_rate_plastic( - const double fluidity, - const double friction) +double +shear_strain_rate_plastic(const double fluidity, const double friction) { return fluidity*friction; } -void compute_shear_strain_rate_plastic(struct simulation* sim) +void +compute_shear_strain_rate_plastic(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) + int i; + for (i=0; i<sim->nz; ++i) sim->gamma_dot_p[i] = shear_strain_rate_plastic(sim->g_ghost[idx1g(i)], sim->mu[i]); } -void compute_shear_velocity(struct simulation* sim) +void +compute_shear_velocity(struct simulation* sim) { + int i; + // TODO: implement iterative solver // Dirichlet BC at bottom sim->v_x[0] = sim->v_x_bot; - for (int i=1; i<sim->nz; ++i) + for (i=1; i<sim->nz; ++i) sim->v_x[i] = sim->v_x[i-1] + sim->gamma_dot_p[i]*sim->dz; } -void compute_effective_stress(struct simulation* sim) +void +compute_effective_stress(struct simulation* sim) { + int i; if (sim->fluid) - for (int i=0; i<sim->nz; ++i) + for (i=0; i<sim->nz; ++i) sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[idx1g(i)]; else - for (int i=0; i<sim->nz; ++i) + for (i=0; i<sim->nz; ++i) sim->sigma_n_eff[i] = sim->sigma_n[i]; } -double cooperativity_length(const double A, - const double d, - const double mu, - const double mu_s) +double +cooperativity_length(const double A, + const double d, + const double mu, + const double mu_s) { return A*d/sqrt(fabs(mu - mu_s)); } -void compute_cooperativity_length(struct simulation* sim) +void +compute_cooperativity_length(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) + int i; + for (i=0; i<sim->nz; ++i) sim->xi[i] = cooperativity_length(sim->A, sim->d, sim->mu[i], sim->mu_s); } -double local_fluidity(const double p, - const double mu, - const double mu_s, - const double b, - const double rho_s, - const double d) +double +local_fluidity(const double p, + const double mu, + const double mu_s, + const double b, + const double rho_s, + const double d) { if (mu <= mu_s) return 0.0; @@ -277,9 +293,11 @@ double local_fluidity(const double p, return sqrt(p/rho_s*d*d) * (mu - mu_s)/(b*mu); } -void compute_local_fluidity(struct simulation* sim) +void +compute_local_fluidity(struct simulation* sim) { - for (int i=0; i<sim->nz; ++i) + int i; + for (i=0; i<sim->nz; ++i) sim->g_ghost[idx1g(i)] = local_fluidity(sim->sigma_n_eff[i], sim->mu[i], sim->mu_s, @@ -288,7 +306,8 @@ void compute_local_fluidity(struct simulation* sim) sim->d); } -void set_bc_neumann(double* g_ghost, const int nz, const int boundary) +void +set_bc_neumann(double* g_ghost, const int nz, const int boundary) { if (boundary == -1) g_ghost[0] = g_ghost[1]; @@ -300,10 +319,11 @@ void set_bc_neumann(double* g_ghost, const int nz, const int boundary) } } -void set_bc_dirichlet(double* g_ghost, - const int nz, - const int boundary, - const double value) +void +set_bc_dirichlet(double* g_ghost, + const int nz, + const int boundary, + const double value) { if (boundary == -1) g_ghost[0] = value; @@ -315,21 +335,25 @@ void set_bc_dirichlet(double* g_ghost, } } -void poisson_solver_1d_cell_update(int i, - const double* g_in, - double* g_out, - double* r_norm, - const double dz, - const double* mu, - const double* p, - const double* xi, - const double mu_s, - const double b, - const double rho_s, - const double d) +void +poisson_solver_1d_cell_update(int i, + const double* g_in, + double* g_out, + double* r_norm, + const double dz, + const double* mu, + const double* p, + const double* xi, + const double mu_s, + const double b, + const double rho_s, + const double d) { - double coorp_term = dz*dz/(2.0*pow(xi[i], 2.0)); - int gi = idx1g(i); + double coorp_term; + int gi; + + coorp_term = dz*dz/(2.0*pow(xi[i], 2.0)); + gi = idx1g(i); g_out[gi] = 1.0/(1.0 + coorp_term)*(coorp_term* local_fluidity(p[i], mu[i], mu_s, b, rho_s, d) + g_in[gi+1]/2.0 @@ -347,15 +371,18 @@ void poisson_solver_1d_cell_update(int i, #endif } -int implicit_1d_jacobian_poisson_solver(struct simulation* sim, - const int max_iter, - const double rel_tol) +int +implicit_1d_jacobian_poisson_solver(struct simulation* sim, + const int max_iter, + const double rel_tol) { - double* g_ghost_out = empty(sim->nz+2); - double* r_norm = empty(sim->nz); + double r_norm_max; + double *g_ghost_out, *r_norm; + int iter, i; - int iter; - double r_norm_max = NAN; + r_norm_max = NAN; + g_ghost_out = empty(sim->nz+2); + r_norm = empty(sim->nz); for (iter=0; iter<max_iter; ++iter) { #ifdef DEBUG @@ -368,7 +395,7 @@ int implicit_1d_jacobian_poisson_solver(struct simulation* sim, /* Neumann BCs resemble free surfaces */ /* set_bc_neumann(sim->g_ghost, sim->nz, +1); */ - for (int i=0; i<sim->nz; ++i) + for (i=0; i<sim->nz; ++i) poisson_solver_1d_cell_update(i, sim->g_ghost, g_ghost_out, @@ -403,23 +430,28 @@ int implicit_1d_jacobian_poisson_solver(struct simulation* sim, return 1; } -void write_output_file(struct simulation* sim, const int normalize) +void +write_output_file(struct simulation* sim, const int normalize) { + int i; char outfile[200]; FILE *fp; + double *v_x_out; + double max_v_x; + sprintf(outfile, "%s.output%05d.txt", sim->name, sim->n_file++); - double *v_x_out = malloc(sim->nz*sizeof(double)); + v_x_out = malloc(sim->nz*sizeof(double)); copy_values(sim->v_x, v_x_out, sim->nz); if (normalize) { - double max_v_x = max(v_x_out, sim->nz); - for (int i=0; i<sim->nz; ++i) + max_v_x = max(v_x_out, sim->nz); + for (i=0; i<sim->nz; ++i) v_x_out[i] /= max_v_x; } fp = fopen(outfile, "w"); if (sim->fluid) - for (int i=0; i<sim->nz; ++i) + for (i=0; i<sim->nz; ++i) fprintf(fp, "%.17g\t%.17g\t%.17g\t%.17g\t%.17g\n", sim->z[i], v_x_out[i],
