commit f463d4d2404fe59625621e49a7db22b70a9202d4
parent 461dca87e55ffa2b46791b589f65cb3b050b533f
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 5 Apr 2019 14:33:35 +0200
Implement implicit Poisson solver
Diffstat:
5 files changed, 233 insertions(+), 29 deletions(-)
diff --git a/arrays.c b/arrays.c
@@ -1,9 +1,10 @@
#include <stdlib.h>
+#include <math.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
+/* 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)
@@ -11,8 +12,8 @@ unsigned int idx3(
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
+/* 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)
@@ -20,22 +21,28 @@ unsigned int idx3g(
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
+/* 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
+/* 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]
+/* Translate a i index in grid with a padding of single into a linear index */
+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 *x = malloc(n*sizeof(double));
@@ -45,7 +52,7 @@ double* linspace(const double lower, const double upper, const int n)
return x;
}
-// Return an array of `n` values with the value 0.0
+/* Return an array of `n` values with the value 0.0 */
double* zeros(const double n)
{
double *x = malloc(n*sizeof(double));
@@ -54,7 +61,7 @@ double* zeros(const double n)
return x;
}
-// Return an array of `n` values with the value 1.0
+/* Return an array of `n` values with the value 1.0 */
double* ones(const double n)
{
double *x = malloc(n*sizeof(double));
@@ -63,3 +70,29 @@ double* ones(const double n)
return x;
}
+/* Return an array of `n` uninitialized values */
+double* empty(const double n)
+{
+ return malloc(n*sizeof(double));
+}
+
+/* Return largest value in array `a` with size `n` */
+double max(const double* a, const int n)
+{
+ double maxval = -INFINITY;
+ for (int 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 minval = +INFINITY;
+ for (int i=0; i<n; ++i)
+ if (a[i] < minval)
+ minval = a[i];
+ return minval;
+}
+
diff --git a/arrays.h b/arrays.h
@@ -15,8 +15,14 @@ unsigned int idx2(
unsigned int idx2g(
const unsigned int i, const unsigned int j, const unsigned int nx);
+unsigned int idx1g(const unsigned int i);
+
double* linspace(const double lower, const double upper, const int n);
double* zeros(const double n);
double* ones(const double n);
+double* empty(const double n);
+
+double max(const double* a, const int n);
+double min(const double* a, const int n);
#endif
diff --git a/damsgaard2013-fig8.png b/damsgaard2013-fig8.png
Binary files differ.
diff --git a/simulation.c b/simulation.c
@@ -1,3 +1,4 @@
+#include <stdio.h>
#include <stdlib.h>
#include "arrays.h"
#include "simulation.h"
@@ -5,11 +6,13 @@
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 */
+ 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->xi = zeros(sim->nz); /* cooperativity length */
+ sim->gamma_dot_p = zeros(sim->nz); /* local shear velocity */
+ 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)
@@ -17,16 +20,160 @@ void free_arrays(struct simulation* sim)
free(sim->z);
free(sim->mu);
free(sim->p);
+ free(sim->xi);
+ free(sim->gamma_dot_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)
+double shear_strain_rate_plastic(
+ const double fluidity,
+ const double friction)
{
- for (int i=0; i<n; ++i)
- shear_strain_rate[i] = fluidity[i]*friction[i];
+ return fluidity*friction;
+}
+
+void compute_shear_strain_rate_plastic(struct simulation* sim)
+{
+ for (int i=0; i<sim->nz; ++i)
+ sim->gamma_dot_p[i] =
+ shear_strain_rate_plastic(sim->g_ghost[idx1g(i)], sim->mu[i]);
+}
+
+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)
+{
+ for (int 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)
+{
+ if (mu <= mu_s)
+ return 0.0;
+ else
+ return sqrt(p/rho_s*d*d) * (mu - mu_s)/(b*mu);
+}
+
+void compute_local_fluidity(struct simulation* sim)
+{
+ for (int i=0; i<sim->nz; ++i)
+ sim->g_ghost[idx1g(i)] = local_fluidity(
+ sim->p[i], sim->mu[i], sim->mu_s, sim->b, sim->rho_s, sim->d);
+}
+
+void set_bc_neumann(double* g_ghost, const int nz, const int boundary)
+{
+ if (boundary == -1)
+ g_ghost[0] = g_ghost[1];
+ else if (boundary == +1)
+ g_ghost[idx1g(nz)+1] = g_ghost[idx1g(nz)];
+ else {
+ fprintf(stderr, "set_bc_neumann: Unknown boundary %d\n", boundary);
+ exit(1);
+ }
+}
+
+void set_bc_dirichlet(
+ double* g_ghost,
+ const int nz,
+ const int boundary,
+ const double value)
+{
+ if (boundary == -1)
+ g_ghost[0] = value;
+ else if (boundary == +1)
+ g_ghost[idx1g(nz)+1] = value;
+ else {
+ fprintf(stderr, "set_bc_dirichlet: Unknown boundary %d\n", boundary);
+ exit(1);
+ }
+}
+
+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));
+ g_out[idx1g(i)] = 1.0/(1.0 + coorp_term)*(coorp_term*
+ local_fluidity(p[i], mu[i], mu_s, b, rho_s, d)
+ + g_in[idx1g(i)+1]/2.0
+ + g_in[idx1g(i)-1]/2.0);
+}
+
+int implicit_1d_jacobian_poisson_solver(
+ struct simulation* sim,
+ const int max_iter,
+ const double rel_tol)
+{
+ double* g_ghost_out = empty(sim->nz);
+ double* r_norm = empty(sim->nz);
+
+ set_bc_neumann(sim->g_ghost, sim->nz, +1);
+ set_bc_dirichlet(sim->g_ghost, sim->nz, -1, 0.0);
+
+ int iter;
+ double r_norm_max = NAN;
+
+ for (iter=0; iter<max_iter; ++iter) {
+
+ for (int i=0; i<sim->nz; ++i)
+ poisson_solver_1d_cell_update(
+ i,
+ sim->g_ghost,
+ g_ghost_out,
+ r_norm,
+ sim->dz,
+ sim->mu,
+ sim->p,
+ sim->xi,
+ sim->mu_s,
+ sim->b,
+ sim->rho_s,
+ sim->d);
+ r_norm_max = max(r_norm, sim->nz);
+
+ double* tmp = sim->g_ghost;
+ sim->g_ghost = g_ghost_out;
+ g_ghost_out = tmp;
+
+ if (r_norm_max <= rel_tol) {
+ free(g_ghost_out);
+ free(r_norm);
+ printf(".. Solution converged after %d iterations\n", iter);
+ return 0;
+ }
+ }
+
+ free(g_ghost_out);
+ free(r_norm);
+ fprintf(stderr, "implicit_1d_jacobian_poisson_solver: ");
+ fprintf(stderr, "Solution did not converge after %d iterations\n", iter);
+ fprintf(stderr, ".. Residual normalized error: %f\n", r_norm_max);
+ return 1;
}
diff --git a/simulation.h b/simulation.h
@@ -40,27 +40,45 @@ struct simulation {
/* nodes along z */
int nz;
- /* origo of axis */
+ /* origo of axis [m] */
double origo_z;
- /* length of domain */
+ /* length of domain [m] */
double L_z;
/* array of cell coordinates */
double* z;
- /* cell spacin */
+ /* cell spacing [m] */
double dz;
/* other arrays */
- double* mu;
- double* p;
- double* v_x;
- double* g_ghost;
+ double* mu; /* static yield friction [-] */
+ double* p; /* effective normal pressure [Pa] */
+ double* xi; /* cooperativity length */
+ double* gamma_dot_p; /* plastic shear strain rate [1/s] */
+ double* v_x; /* shear velocity [m/s] */
+ double* g_ghost; /* fluidity with ghost nodes */
};
void prepare_arrays(struct simulation* sim);
void free_arrays(struct simulation* sim);
+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);
+
+int implicit_1d_jacobian_poisson_solver(
+ struct simulation* sim,
+ const int max_iter,
+ const double rel_tol);
+
#endif
