tadd untested grain, array, and util files - granular - granular dynamics simulation
HTML git clone git://src.adamsgaard.dk/granular
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
DIR commit 4a58873534515272d0df10f6f023727a554bce11
DIR parent 0a66564a557db790b602d756741d798909a76f71
HTML Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 18 Mar 2021 00:00:54 +0100
add untested grain, array, and util files
Diffstat:
A arrays.c | 292 +++++++++++++++++++++++++++++++
A arrays.h | 54 +++++++++++++++++++++++++++++++
A grain.c | 147 +++++++++++++++++++++++++++++++
A grain.h | 39 +++++++++++++++++++++++++++++++
A simulation.h | 162 ++++++++++++++++++++++++++++++
A util.c | 78 +++++++++++++++++++++++++++++++
A util.h | 10 ++++++++++
7 files changed, 782 insertions(+), 0 deletions(-)
---
DIR diff --git a/arrays.c b/arrays.c
t@@ -0,0 +1,292 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define DEG2RAD(x) (x*M_PI/180.0)
+
+void
+check_magnitude(const char *func_name, int limit, int value)
+{
+ if (value < limit) {
+ fprintf(stderr, "error: %s: input size %d is less than %d\n",
+ func_name, value, limit);
+ exit(1);
+ }
+}
+
+/* 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);
+}
+
+/* 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)
+{
+ int i;
+ double *x;
+ double dx;
+
+ check_magnitude(__func__, 1, n);
+ 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;
+
+ check_magnitude(__func__, 2, n);
+ 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)
+{
+ int i;
+ double *x;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+ double *x;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+ double *x;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ check_magnitude(__func__, 1, n);
+ return malloc(n * sizeof(double));
+}
+
+/* Return largest value in array `a` with size `n` */
+double
+max(const double *a, const int n)
+{
+ int i;
+ double maxval;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+ double minval;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+
+ check_magnitude(__func__, 1, n);
+ for (i = 0; i < n; ++i)
+ printf("%.17g\n", a[i]);
+}
+
+void
+print_arrays(const double *a, const double *b, const int n)
+{
+ int i;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+ double max_b;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+
+ check_magnitude(__func__, 1, n);
+ 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)
+{
+ int i;
+
+ check_magnitude(__func__, 1, n);
+ for (i = 0; i < n; ++i)
+ out[i] = in[i];
+}
+
+double *
+copy(const double *in, const int n)
+{
+ double *out;
+
+ check_magnitude(__func__, 1, n);
+ out = empty(n);
+ copy_values(in, out, n);
+ return out;
+}
+
+double *
+normalize(const double *in, const int n)
+{
+ int i;
+ double max_val;
+ double *out;
+
+ check_magnitude(__func__, 1, n);
+ out = malloc(n * sizeof(double));
+ copy_values(in, out, n);
+ max_val = max(out, n);
+
+ if (max_val == 0.0)
+ max_val = 1.0;
+
+ for (i = 0; i < n; ++i)
+ out[i] /= max_val;
+
+ return out;
+}
DIR diff --git a/arrays.h b/arrays.h
t@@ -0,0 +1,54 @@
+#include <stdio.h>
+
+#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);
+
+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);
+double * initval(const double value, const int n);
+double * empty(const int n);
+
+double max(const double *a, const int n);
+double min(const double *a, const int n);
+
+void print_array(const double *a, const int n);
+void print_arrays(const double *a, const double *b, const int n);
+void print_arrays_2nd_normalized(const double *a, const double *b, const int n);
+void print_three_arrays(
+ const double *a,
+ const double *b,
+ const double *c,
+ const int n);
+
+void fprint_arrays(FILE *fp, const double *a, const double *b, const int n);
+
+void fprint_three_arrays(
+ FILE *fp,
+ const double *a,
+ const double *b,
+ const double *c,
+ const int n);
+
+void copy_values(const double *in, double *out, const int n);
+double * copy(const double *in, const int n);
+double * normalize(const double *in, const int n);
+
+#endif
DIR diff --git a/grain.c b/grain.c
t@@ -0,0 +1,147 @@
+#include <stdio.h>
+#include "util.h"
+#include "grain.h"
+
+extern size_t nd;
+
+void
+grain_defaults(struct grain *g)
+{
+ size_t i;
+
+ g->radius = 1.0;
+ for (i = 0; i < nd; i++) {
+ g->pos[i]
+ = g->vel[i]
+ = g->acc[i]
+ = g->force[i]
+ = g->angpos[i]
+ = g->angvel[i]
+ = g->angacc[i]
+ = g->torque[i]
+ = g->disp[i]
+ = g->forceext[i]
+ = g->contact_stress[i]
+ = 0.0;
+ g->gridpos[i] = 0;
+ }
+ g->fixed = 0;
+ g->rotating = 1;
+ g->enabled = 1;
+ g->youngs_modulus = 1e9;
+ g->poissons_ratio = 0.2;
+ g->friction_coeff = 0.4;
+ g->tensile_strength = 1e8;
+ g->ncontacts = 0;
+ g->thermal_energy = 0.0;
+ g->color = 0;
+}
+
+static void
+print_padded_nd_double(FILE *stream, double *arr)
+{
+ size_t i;
+
+ for (i = 0; i < 3; i++)
+ if (i < nd)
+ fprintf(stream, "%.17g\t", arr[i]);
+ else
+ fprintf(stream, "0.0\t");
+}
+
+static void
+print_padded_nd_int(FILE *stream, size_t *arr)
+{
+ size_t i;
+
+ for (i = 0; i < 3; i++)
+ if (i < nd)
+ fprintf(stream, "%d\t", arr[i]);
+ else
+ fprintf(stream, "0.0\t");
+}
+
+void
+print_grain(FILE *stream, const struct grain *g)
+{
+ size_t i;
+
+ fprintf(stream, "%d\t", g->radius);
+ print_padded_nd_double(stream, g->pos);
+ print_padded_nd_double(stream, g->vel);
+ print_padded_nd_double(stream, g->acc);
+ print_padded_nd_double(stream, g->force);
+ print_padded_nd_double(stream, g->angpos);
+ print_padded_nd_double(stream, g->angvel);
+ print_padded_nd_double(stream, g->angacc);
+ print_padded_nd_double(stream, g->torque);
+ print_padded_nd_double(stream, g->disp);
+ print_padded_nd_double(stream, g->forceext);
+ fprintf(stream, "%.17g\t", g->density);
+ fprintf(stream, "%d\t", g->fixed);
+ fprintf(stream, "%d\t", g->rotating);
+ fprintf(stream, "%d\t", g->enabled);
+ fprintf(stream, "%.17g\t", g->youngs_modulus);
+ fprintf(stream, "%.17g\t", g->poissons_ratio);
+ fprintf(stream, "%.17g\t", g->friction_coefficient);
+ fprintf(stream, "%.17g\t", g->tensile_strength);
+ fprintf(stream, "%.17g\t", g->shear_strength);
+ fprintf(stream, "%.17g\t", g->fracture_toughness);
+ print_padded_nd_int(stream, g->grid_position);
+ fprintf(stream, "%d\t", g->ncontacts);
+ fprintf(stream, "%.17g\t", g->granular_stress);
+ fprintf(stream, "%.17g\t", g->thermal_energy);
+ fprintf(stream, "%d\n", g->color);
+}
+
+int
+check_grain_values(const struct grain *g)
+{
+ int i;
+ int status = 0;
+
+ check_float_positive("grain->radius", g->radius, &status);
+
+ for (i = 0; i < nd; i++) {
+ check_float("grain->pos", g->pos[i], &status);
+ check_float("grain->vel", g->vel[i], &status);
+ check_float("grain->acc", g->acc[i], &status);
+ check_float("grain->force", g->force[i], &status);
+ check_float("grain->angpos", g->angpos[i], &status);
+ check_float("grain->angvel", g->angvel[i], &status);
+ check_float("grain->angacc", g->angacc[i], &status);
+ check_float("grain->torque", g->torque[i], &status);
+ check_float("grain->disp", g->disp[i], &status);
+ check_float("grain->forceext", g->forceext[i], &status);
+ check_float("grain->contact_stress", g->contact_stress[i], &status);
+ }
+
+ check_bool("grain->fixed", g->fixed, &status);
+ check_bool("grain->rotating", g->rotating, &status);
+ check_bool("grain->enabled", g->enabled, &status);
+
+ check_float_non_negative("grain->youngs_modulus",
+ g->youngs_modulus,
+ &status);
+ check_float_non_negative("grain->poissons_ratio",
+ g->youngs_modulus,
+ &status);
+ check_float_non_negative("grain->friction_coeff",
+ g->youngs_modulus,
+ &status);
+ check_float_non_negative("grain->tensile_strength",
+ g->youngs_modulus,
+ &status);
+
+ for (i = 0; i < nd; i++)
+ if (g->gridpos[i] < 0)
+ warn_parameter_value("grain->gridpos is not 0 or 1",
+ (double)g->gridpos[i], &status);
+
+ if (g->ncontacts[i] < 0)
+ warn_parameter_value("grain->ncontacts is not 0 or 1",
+ (double)g->ncontacts[i], &status);
+ check_float_non_negative("grain->thermal_energy", g->thermal_energy, &status);
+
+ return status;
+}
DIR diff --git a/grain.h b/grain.h
t@@ -0,0 +1,39 @@
+#ifndef GRANULAR_GRAIN_
+#define GRANULAR_GRAIN_
+
+extern size_t nd;
+
+struct grain {
+
+ double radius;
+ double pos[nd];
+ double vel[nd];
+ double acc[nd];
+ double force[nd];
+ double angpos[nd];
+ double angvel[nd];
+ double angacc[nd];
+ double torque[nd];
+ double disp[nd];
+ double forceext[nd];
+ double density;
+ size_t fixed;
+ size_t rotating;
+ size_t enabled;
+ double youngs_modulus;
+ double poissons_ratio;
+ double friction_coeff;
+ double tensile_strength;
+ size_t gridpos[nd];
+ size_t ncontacts;
+ double contact_stress[nd];
+ double thermal_energy;
+ size_t color;
+};
+
+void init_grain(struct grain *g);
+void print_grain(FILE *stream, const struct grain *g);
+int check_grain_values(const struct grain *g);
+double grain_volume(const struct grain *g);
+
+#endif
DIR diff --git a/simulation.h b/simulation.h
t@@ -0,0 +1,162 @@
+#ifndef SIMULATION_
+#define SIMULATION_
+
+#include "arrays.h"
+
+#define DEFAULT_SIMULATION_NAME "unnamed_simulation"
+#define PI 3.14159265358979323846
+#define DEG2RAD(x) (x*PI/180.0)
+
+extern struct simulation sim;
+
+/* Simulation settings */
+struct simulation {
+
+ /* simulation name to use for output files */
+ char name[100];
+
+ /* gravitational acceleration magnitude [m/s^2] */
+ double G;
+
+ /* normal stress from the top wall [Pa] */
+ double P_wall;
+
+ /* optionally fix top shear velocity to this value [m/s] */
+ double v_x_fix;
+
+ /* optionally fix top shear velocity to this value [m/s] */
+ double v_x_limit;
+
+ /* 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;
+
+ /* material cohesion [Pa] */
+ double C;
+
+ /* representative grain size [m] */
+ double d;
+
+ /* grain material density [kg/m^3] */
+ double rho_s;
+
+ /* nodes along z */
+ int nz;
+
+ /* origo of axis [m] */
+ double origo_z;
+
+ /* length of domain [m] */
+ double L_z;
+
+ /* array of cell coordinates */
+ double *z;
+
+ /* cell spacing [m] */
+ double dz;
+
+ /* current time [s] */
+ double t;
+
+ /* end time [s] */
+ double t_end;
+
+ /* time step length [s] */
+ double dt;
+
+ /* interval between output files [s] */
+ double file_dt;
+
+ /* output file number */
+ int n_file;
+
+ double transient;
+ double phi_min;
+ double phi_max;
+ double dilatancy_constant;
+
+ /* Fluid parameters */
+ int fluid; /* flag to switch fluid on (1) or off (0) */
+ double p_f_top; /* fluid pressure at the top [Pa] */
+ double p_f_mod_ampl; /* amplitude of fluid pressure variations [Pa] */
+ double p_f_mod_freq; /* frequency of fluid pressure variations [s^-1] */
+ double p_f_mod_phase; /* phase of fluid pressure variations [s^-1] */
+ double p_f_mod_pulse_time; /* single pressure pulse at this time [s] */
+ int p_f_mod_pulse_shape; /* waveform for fluid-pressure pulse */
+ double beta_f; /* adiabatic fluid compressibility [Pa^-1] */
+ double alpha; /* adiabatic grain compressibility [Pa^-1] */
+ double mu_f; /* fluid dynamic viscosity [Pa*s] */
+ double rho_f; /* fluid density [kg/m^3] */
+ double D; /* diffusivity [m^2/s], overrides k, beta_f, alpha, mu_f */
+
+ /* arrays */
+ double *mu; /* static yield friction [-] */
+ double *mu_c; /* critical-state static yield friction [-] */
+ double *sigma_n_eff; /* effective normal pressure [Pa] */
+ double *sigma_n; /* normal stress [Pa] */
+ double *p_f_ghost; /* fluid pressure [Pa] */
+ double *p_f_dot; /* fluid pressure change [Pa/s] */
+ double *p_f_dot_expl; /* fluid pressure change (explicit solution) [Pa/s] */
+ double *p_f_dot_impl; /* fluid pressure change (implicit solution) [Pa/s] */
+ double *p_f_dot_impl_r_norm; /* normalized residual fluid pressure change [-] */
+ double *k; /* hydraulic permeability [m^2] */
+ double *phi; /* porosity [-] */
+ double *phi_c; /* critical-state porosity [-] */
+ double *phi_dot; /* porosity change [s^-1] */
+ double *xi; /* cooperativity length */
+ double *gamma_dot_p; /* plastic shear strain rate [s^-1] */
+ double *v_x; /* shear velocity [m/s] */
+ double *g_local; /* local fluidity */
+ double *g_ghost; /* fluidity with ghost nodes */
+ double *g_r_norm; /* normalized residual of fluidity field */
+ double *I; /* inertia number [-] */
+ double *tan_psi; /* tan(dilatancy_angle) [-] */
+ double *old_val; /* temporary storage for iterative solvers */
+ double *fluid_old_val;/* temporary storage for fluid iterative solver */
+ double *tmp_ghost; /* temporary storage for iterative solvers */
+ double *p_f_dot_old; /* temporary storage for old p_f_dot */
+};
+
+void init_sim(struct simulation *sim);
+void prepare_arrays(struct simulation *sim);
+void free_arrays(struct simulation *sim);
+void check_simulation_parameters(struct simulation *sim);
+void lithostatic_pressure_distribution(struct simulation *sim);
+void compute_effective_stress(struct simulation *sim);
+
+void set_bc_neumann(double *a,
+ const int nz,
+ const int boundary,
+ const double df,
+ const double dx);
+
+void set_bc_dirichlet(double *a,
+ const int nz,
+ const int boundary,
+ const double value);
+
+double residual(double new_val, double old_val);
+
+void write_output_file(struct simulation *sim, const int normalize);
+void print_output(struct simulation *sim, FILE *fp, const int normalize);
+
+int coupled_shear_solver(struct simulation *sim,
+ const int max_iter,
+ const double rel_tol);
+
+void set_coupled_fluid_transient_timestep(struct simulation *sim, const double safety);
+
+double find_flux(const struct simulation *sim);
+
+#endif
DIR diff --git a/util.c b/util.c
t@@ -0,0 +1,78 @@
+#include <stdio.h>
+
+void
+warn_parameter_value(const char message[],
+ const double value,
+ int *status)
+{
+ fprintf(stderr,
+ "%s: %s (%.17g)\n",
+ __func__, message, value);
+ *status = 1;
+}
+
+void
+check_float(const char name[], const double value, int *status)
+{
+ if (isnan(value)) {
+ char message[100];
+
+ snprintf(message, sizeof(message), "%s is NaN", name);
+ warn_parameter_value(message, value, status);
+ *status = 1;
+ } else if (isinf(value)) {
+ char message[100];
+
+ snprintf(message, sizeof(message), "%s is infinite", name);
+ warn_parameter_value(message, value, status);
+ *status = 1;
+ }
+}
+
+void
+check_float_non_negative(const char name[], const double value, int *status)
+{
+ char message[100];
+
+ check_float(*name, value, *status);
+ if (value < 0.0) {
+ snprintf(message, sizeof(message), "%s is negative", name);
+ warn_parameter_value(message, value, status);
+ *status = 1;
+ }
+}
+
+void
+check_float_positive(const char name[], const double value, int *status)
+{
+ char message[100];
+
+ check_float(*name, value, *status);
+ if (value <= 0.0) {
+ snprintf(message, sizeof(message), "%s is not positive", name);
+ warn_parameter_value(message, value, status);
+ *status = 1;
+ }
+}
+
+void
+check_bool(const char name[], const int value, int *status)
+{
+ if (value < 0 || value > 1)
+ warn_parameter_value("%s is not 0 or 1",
+ (double)value,
+ &status);
+ *status = 1;
+}
+
+void
+check_int_non_negative(const char name[], const int value, int *status)
+{
+ char message[100];
+
+ if (value < 0) {
+ snprintf(message, sizeof(message), "%s is negative", name);
+ warn_parameter_value(message, (double)value, status);
+ *status = 1;
+ }
+}
DIR diff --git a/util.h b/util.h
t@@ -0,0 +1,10 @@
+#ifndef GRANULAR_UTIL_
+#define GRANULAR_UTIL_
+
+void warn_parameter_value(const char message[], const double value, int *status);
+void check_float(const char name[], const double value, int *status);
+void check_float_non_negative(const char name[], const double value, int *status);
+void check_float_positive(const char name[], const double value, int *status);
+void check_bool(const char name[], const double value, int *status);
+void check_int_non_negative(const char name[], const int value, int *status);
+#endif