1d_fd_simple_shear

commit bf84e746592960c537206a8913ffbef07110f228
parent 16ea3677eab3418adb201a5c2e88cac7761eded5
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Fri,  5 Apr 2019 15:19:28 +0200

Add debug information and compute missing values

Diffstat:
M
1d_fd_simple_shear_damsgaard2013.h
 | 
2
+-
M
Makefile
 | 
4
++++
M
arrays.c
 | 
6
++++++
M
arrays.h
 | 
2
++
M
main.c
 | 
18
++++++++++++++++++
M
simulation.c
 | 
16
+++++++++++++---
M
simulation.h
 | 
3
+++

7 files changed, 47 insertions(+), 4 deletions(-)
diff --git a/1d_fd_simple_shear_damsgaard2013.h b/1d_fd_simple_shear_damsgaard2013.h
@@ -18,7 +18,7 @@ struct simulation init_sim(void)
     sim.P_wall = 10e3; /* larger normal stress deepens the shear depth */
     sim.mu_wall = 0.40;
 
-    sim.nz = 100;
+    sim.nz = 10;
 
     /* lower values of A mean that the velocity curve can have sharper curves,
      * e.g. at the transition from μ ≈ μ_s */
diff --git a/Makefile b/Makefile
@@ -17,6 +17,10 @@ default: 1d_fd_simple_shear
 1d_fd_simple_shear: $(OBJ) $(HDR)
 	$(CC) $(LDFLAGS) $(OBJ) -o $@
 
+.PHONY: watch
+watch:
+	echo $(SRC) $(HDR) | tr ' ' '\n' | entr -s 'make && ./1d_fd_simple_shear'
+
 clean:
 	$(RM) *.png
 	$(RM) *.o
diff --git a/arrays.c b/arrays.c
@@ -1,4 +1,5 @@
 #include <stdlib.h>
+#include <stdio.h>
 #include <math.h>
 
 #define DEG2RAD(x) (x*M_PI/180.0)
@@ -96,3 +97,8 @@ double min(const double* a, const int n)
     return minval;
 }
 
+void print_array(const double* a, const int n)
+{
+    for (int i=0; i<n; ++i)
+        printf("%g\n", a[i]);
+}
diff --git a/arrays.h b/arrays.h
@@ -25,4 +25,6 @@ double* empty(const double 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);
+
 #endif
diff --git a/main.c b/main.c
@@ -1,4 +1,5 @@
 #include <stdio.h>
+#include <stdlib.h>
 #include <math.h>
 
 #include "simulation.h"
@@ -13,6 +14,23 @@ int main(int argc, char** argv) {
     struct simulation sim = init_sim();
     prepare_arrays(&sim);
 
+    init_pressure(&sim);
+    init_friction(&sim);
+    compute_cooperativity_length(&sim);
+
+    puts("\n## Before solver");
+    puts(".. p:"); print_array(sim.p, sim.nz);
+    puts(".. mu:"); print_array(sim.mu, sim.nz);
+
+    if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5))
+        exit(1);
+
+    compute_shear_strain_rate_plastic(&sim);
+
+    puts("\n## After solver");
+    puts(".. g:"); print_array(sim.g_ghost, sim.nz+2);
+    puts(".. v_x:"); print_array(sim.v_x, sim.nz);
+
     free_arrays(&sim);
     return 0;
 }
diff --git a/simulation.c b/simulation.c
@@ -120,10 +120,20 @@ void poisson_solver_1d_cell_update(
         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*
+    int 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[idx1g(i)+1]/2.0
-            + g_in[idx1g(i)-1]/2.0);
+            + g_in[gi+1]/2.0
+            + g_in[gi-1]/2.0);
+
+    r_norm[i] = pow(g_out[gi] - g_in[gi], 1.0) / (pow(g_out[gi], 2.0) + 1e-16);
+
+    printf("-- %d --------------\n", i);
+    printf("coorp_term: %g\n", coorp_term);
+    printf("   g_local: %g\n", local_fluidity(p[i], mu[i], mu_s, b, rho_s, d));
+    printf("      g_in: %g\n", g_in[gi]);
+    printf("     g_out: %g\n", g_out[gi]);
+    printf("    r_norm: %g\n", r_norm[gi]);
 }
 
 int implicit_1d_jacobian_poisson_solver(
diff --git a/simulation.h b/simulation.h
@@ -76,6 +76,9 @@ void set_bc_dirichlet(
         const int boundary,
         const double value);
 
+void compute_cooperativity_length(struct simulation* sim);
+void compute_shear_strain_rate_plastic(struct simulation* sim);
+
 int implicit_1d_jacobian_poisson_solver(
         struct simulation* sim,
         const int max_iter,