tMerge branch 'debug_diffusion' into 'master' - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
HTML git clone git://src.adamsgaard.dk/cngf-pf
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
DIR commit edf97b34d7ba2dcd925192cd64479d16c4767b2a
DIR parent b2f8d047df70db4aa5593bfafe9bd594312e1481
HTML 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:
M .gitlab-ci.yml | 10 ++++------
D 1d_fd_simple_shear_henann_kamrin20… | 36 -------------------------------
M Makefile | 113 +++----------------------------
M README.md | 10 +++++-----
M arrays.c | 174 ++++++++++++++++++++-----------
M arrays.h | 1 +
D diurnal.gif | 0
R damsgaard2013-fig8.png -> doc/dams… | 0
R iverson2010-fig2a.png -> doc/ivers… | 0
R iverson2010-fig2b.png -> doc/ivers… | 0
R kamb1991-fig1.png -> doc/kamb1991-… | 0
R tulaczyk2006-fig1.png -> doc/tulac… | 0
R 1d_fd_simple_shear.gp -> examples/… | 0
R 1d_fd_simple_shear.png -> examples… | 0
R 1d_fd_simple_shear_fluid.gp -> exa… | 0
R 1d_fd_simple_shear_rheology.gp -> … | 0
R 1d_fd_simple_shear_rheology.png ->… | 0
R 1d_fd_simple_shear_rheology_iverso… | 0
R 1d_fd_simple_shear_rheology_iverso… | 0
R 1d_fd_simple_shear_rheology_kamb.g… | 0
R 1d_fd_simple_shear_rheology_kamb.p… | 0
R 1d_fd_simple_shear_rheology_tulacz… | 0
R 1d_fd_simple_shear_rheology_tulacz… | 0
R BlueSeq.plt -> examples/BlueSeq.plt | 0
A examples/Makefile | 105 +++++++++++++++++++++++++++++++
A examples/diurnal.gif | 0
R diurnal.timeseries.gp -> examples/… | 0
M fluid.c | 173 +++++++++++++++++++------------
M fluid.h | 2 ++
M main.c | 83 +++++++++++++------------------
M simulation.c | 228 ++++++++++++++++++-------------
31 files changed, 516 insertions(+), 419 deletions(-)
---
DIR diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
t@@ -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
DIR diff --git a/1d_fd_simple_shear_henann_kamrin2016.h b/1d_fd_simple_shear_henann_kamrin2016.h
t@@ -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
DIR diff --git a/Makefile b/Makefile
t@@ -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
t@@ -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
DIR diff --git a/README.md b/README.md
t@@ -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
t@@ -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.
-
+
DIR diff --git a/arrays.c b/arrays.c
t@@ -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];
}
DIR diff --git a/arrays.h b/arrays.h
t@@ -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);
DIR diff --git a/diurnal.gif b/diurnal.gif
Binary files differ.
DIR diff --git a/damsgaard2013-fig8.png b/doc/damsgaard2013-fig8.png
Binary files differ.
DIR diff --git a/iverson2010-fig2a.png b/doc/iverson2010-fig2a.png
Binary files differ.
DIR diff --git a/iverson2010-fig2b.png b/doc/iverson2010-fig2b.png
Binary files differ.
DIR diff --git a/kamb1991-fig1.png b/doc/kamb1991-fig1.png
Binary files differ.
DIR diff --git a/tulaczyk2006-fig1.png b/doc/tulaczyk2006-fig1.png
Binary files differ.
DIR diff --git a/1d_fd_simple_shear.gp b/examples/1d_fd_simple_shear.gp
DIR diff --git a/1d_fd_simple_shear.png b/examples/1d_fd_simple_shear.png
Binary files differ.
DIR diff --git a/1d_fd_simple_shear_fluid.gp b/examples/1d_fd_simple_shear_fluid.gp
DIR diff --git a/1d_fd_simple_shear_rheology.gp b/examples/1d_fd_simple_shear_rheology.gp
DIR diff --git a/1d_fd_simple_shear_rheology.png b/examples/1d_fd_simple_shear_rheology.png
Binary files differ.
DIR diff --git a/1d_fd_simple_shear_rheology_iverson.gp b/examples/1d_fd_simple_shear_rheology_iverson.gp
DIR diff --git a/1d_fd_simple_shear_rheology_iverson.png b/examples/1d_fd_simple_shear_rheology_iverson.png
Binary files differ.
DIR diff --git a/1d_fd_simple_shear_rheology_kamb.gp b/examples/1d_fd_simple_shear_rheology_kamb.gp
DIR diff --git a/1d_fd_simple_shear_rheology_kamb.png b/examples/1d_fd_simple_shear_rheology_kamb.png
Binary files differ.
DIR diff --git a/1d_fd_simple_shear_rheology_tulaczyk.gp b/examples/1d_fd_simple_shear_rheology_tulaczyk.gp
DIR diff --git a/1d_fd_simple_shear_rheology_tulaczyk.png b/examples/1d_fd_simple_shear_rheology_tulaczyk.png
Binary files differ.
DIR diff --git a/BlueSeq.plt b/examples/BlueSeq.plt
DIR diff --git a/examples/Makefile b/examples/Makefile
t@@ -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
DIR diff --git a/examples/diurnal.gif b/examples/diurnal.gif
Binary files differ.
DIR diff --git a/diurnal.timeseries.gp b/examples/diurnal.timeseries.gp
DIR diff --git a/fluid.c b/fluid.c
t@@ -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
t@@ -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);
t@@ -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,
t@@ -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]);
t@@ -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);
t@@ -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);
DIR diff --git a/fluid.h b/fluid.h
t@@ -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,
DIR diff --git a/main.c b/main.c
t@@ -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"
t@@ -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"
t@@ -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,
t@@ -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"
t@@ -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'},
t@@ -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 */
t@@ -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",
t@@ -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");
t@@ -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);
t@@ -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],
DIR diff --git a/simulation.c b/simulation.c
t@@ -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,
t@@ -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);
t@@ -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);
t@@ -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)
t@@ -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) {
t@@ -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;
t@@ -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,
t@@ -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];
t@@ -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;
t@@ -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
t@@ -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
t@@ -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,
t@@ -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],