commit 461dca87e55ffa2b46791b589f65cb3b050b533f
parent 4e12789481a70d5f222312b833524be5d19963e2
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 5 Apr 2019 13:32:41 +0200
Add unfinished implementation in C
Diffstat:
7 files changed, 307 insertions(+), 0 deletions(-)
diff --git a/1d_fd_simple_shear_damsgaard2013.h b/1d_fd_simple_shear_damsgaard2013.h
@@ -0,0 +1,69 @@
+#ifndef ONED_FD_SIMPLE_SHEAR_
+#define ONED_FD_SIMPLE_SHEAR_
+
+#include <math.h>
+#include "arrays.h"
+#include "simulation.h"
+
+#define PI 3.14159265358979323846
+#define DEG2RAD(x) (x*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.nz = 100;
+
+ /* lower values of A mean that the velocity curve can have sharper curves,
+ * e.g. at the transition from μ ≈ μ_s */
+ //sim.A = 0.48; /* Henann and Kamrin 2016 */
+ sim.A = 0.40; /* Loose fit to Damsgaard et al 2013 */
+
+ /* lower values of b mean larger shear velocity for a given stress ratio above */
+ /* yield */
+ sim.b = 0.9377; /* Henann and Kamrin 2016 */
+
+ /*sim.mu_s = 0.3819; // Henann and Kamrin 2016 */
+ sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */
+
+ /*sim.phi = 0.38; // Henann and Kamrin 2016 */
+ sim.phi = 0.25; /* Damsgaard et al 2013 */
+
+ /* lower values of d mean that the shear velocity curve can have sharper
+ * curves, e.g. at the transition from μ ≈ μ_s */
+ /*sim.d = 1e-3; // Henann and Kamrin 2016 */
+ sim.d = 0.04; /* Damsgaard et al 2013 */
+
+ /* grain material density [kg/m^3] */
+ /*sim.rho_s = 2.485e3; // Henann and Kamrin 2016 */
+ sim.rho_s = 2.6e3; /* Damsgaard et al 2013 */
+
+ /* Spatial settings */
+ sim.origo_z = 0.0;
+ /*sim.L_z = 20.0*d; // Henann and Kamrin 2016 */
+ sim.L_z = 0.7; /* Damsgaard et al 2013 */
+
+ return sim;
+}
+
+void init_pressure(struct simulation* sim)
+{
+ for (int i=0; i<sim->nz; ++i)
+ sim->p[i] = sim->P_wall +
+ (1.0 - sim->phi)*sim->rho_s*sim->G*(sim->L_z - sim->z[i]);
+}
+
+void init_friction(struct simulation* sim)
+{
+ for (int i=0; i<sim->nz; ++i)
+ sim->mu[i] = sim->mu_wall /
+ (1.0 + (1.0 - sim->phi)*sim->rho_s*sim->G*(sim->L_z - sim->z[i]));
+}
+
+#endif
diff --git a/1d_fd_simple_shear_henann_kamrin2016.h b/1d_fd_simple_shear_henann_kamrin2016.h
@@ -0,0 +1,35 @@
+#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.nz = 100;
+
+ sim.A = 0.48;
+ sim.b = 0.9377;
+ sim.mu_s = 0.3819;
+ sim.phi = 0.38;
+ 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/arrays.c b/arrays.c
@@ -0,0 +1,65 @@
+#include <stdlib.h>
+
+#define DEG2RAD(x) (x*M_PI/180.0)
+
+// 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)
+{
+ 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)
+{
+ 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)
+{
+ 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)
+{
+ return i+1 + (nx+2)*(j+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 *x = malloc(n*sizeof(double));
+ double dx = (upper - lower)/(double)n;
+ for (int i=0; i<n; ++i)
+ x[i] = lower + dx*i;
+ return x;
+}
+
+// Return an array of `n` values with the value 0.0
+double* zeros(const double n)
+{
+ double *x = malloc(n*sizeof(double));
+ for (int 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 double n)
+{
+ double *x = malloc(n*sizeof(double));
+ for (int i=0; i<n; ++i)
+ x[i] = 1.0;
+ return x;
+}
+
diff --git a/arrays.h b/arrays.h
@@ -0,0 +1,22 @@
+#ifndef ARRAYS_
+#define ARRAYS_
+
+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 idx3g(
+ const unsigned int i, const unsigned int j, const unsigned int k,
+ const unsigned int nx, const unsigned int ny);
+
+unsigned int idx2(
+ 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);
+
+double* linspace(const double lower, const double upper, const int n);
+double* zeros(const double n);
+double* ones(const double n);
+
+#endif
diff --git a/main.c b/main.c
@@ -0,0 +1,18 @@
+#include <stdio.h>
+#include <math.h>
+
+#include "simulation.h"
+
+/* set up parameter values for this simulation */
+#include "1d_fd_simple_shear_damsgaard2013.h"
+/* #include "1d_fd_simple_shear_henann_kamrin2016.h" */
+
+
+int main(int argc, char** argv) {
+
+ struct simulation sim = init_sim();
+ prepare_arrays(&sim);
+
+ free_arrays(&sim);
+ return 0;
+}
diff --git a/simulation.c b/simulation.c
@@ -0,0 +1,32 @@
+#include <stdlib.h>
+#include "arrays.h"
+#include "simulation.h"
+
+void prepare_arrays(struct simulation* sim)
+{
+ sim->z = linspace(sim->origo_z, sim->L_z, sim->nz); /* spatial coordinates */
+ sim->dz = sim->z[1] - sim->z[0]; /* cell spacing */
+ sim->mu = zeros(sim->nz); /* local stress ratio */
+ sim->p = zeros(sim->nz); /* local pressure */
+ sim->v_x = zeros(sim->nz); /* local shear velocity */
+ sim->g_ghost = zeros(sim->nz+2); /* local fluidity with ghost nodes */
+}
+
+void free_arrays(struct simulation* sim)
+{
+ free(sim->z);
+ free(sim->mu);
+ free(sim->p);
+ free(sim->v_x);
+ free(sim->g_ghost);
+}
+
+void shear_strain_rate_plastic(
+ const double* fluidity,
+ const double* friction,
+ double* shear_strain_rate,
+ const int n)
+{
+ for (int i=0; i<n; ++i)
+ shear_strain_rate[i] = fluidity[i]*friction[i];
+}
diff --git a/simulation.h b/simulation.h
@@ -0,0 +1,66 @@
+#ifndef SIMULATION_
+#define SIMULATION_
+
+#include <math.h>
+#include "arrays.h"
+
+/* Simulation settings */
+struct simulation {
+
+ /* gravitational acceleration magnitude [m/s^2] */
+ double G;
+
+ /* wall parameters */
+ double P_wall; // normal stress from top wall [Pa]
+
+ /* bottom velocity along x [m/s] */
+ double v_x_bot;
+
+ /* stress ratio at top wall */
+ double mu_wall;
+
+ /* nonlocal amplitude [-] */
+ double A;
+
+ /* rate dependence beyond yield [-] */
+ double b;
+
+ /* bulk and critical state static yield friction coefficient [-] */
+ double mu_s;
+
+ /* porosity [-] */
+ double phi;
+
+ /* representative grain size [m] */
+ double d;
+
+ /* grain material density [kg/m^3] */
+ double rho_s;
+
+ /* nodes along z */
+ int nz;
+
+ /* origo of axis */
+ double origo_z;
+
+ /* length of domain */
+ double L_z;
+
+ /* array of cell coordinates */
+ double* z;
+
+ /* cell spacin */
+ double dz;
+
+ /* other arrays */
+ double* mu;
+ double* p;
+ double* v_x;
+ double* g_ghost;
+};
+
+
+void prepare_arrays(struct simulation* sim);
+void free_arrays(struct simulation* sim);
+
+#endif
