tFirst attempt at dampening the coupled solver. - cngf-pf - continuum model for granular flows with pore-pressure dynamics (renamed from 1d_fd_simple_shear)
 (HTM) git clone git://src.adamsgaard.dk/cngf-pf
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit cd1860a031d74199c9308ee12bee5dea8572f7be
 (DIR) parent 2379b7b9d10d2e3f2e70c3e8ab39c22f2ca88189
 (HTM) Author: Ian Madden <iamadden@stanford.edu>
       Date:   Mon,  8 Mar 2021 20:25:39 -0800
       
       First attempt at dampening the coupled solver.
       
       Signed-off-by: Anders Damsgaard <anders@adamsgaard.dk>
       
       Diffstat:
         M fluid.c                             |       9 +++++++--
         M simulation.c                        |       2 ++
         M simulation.h                        |       1 +
       
       3 files changed, 10 insertions(+), 2 deletions(-)
       ---
 (DIR) diff --git a/fluid.c b/fluid.c
       t@@ -170,13 +170,16 @@ darcy_solver_1d(struct simulation *sim,
                        const double rel_tol)
        {
                int i, iter, solved = 0;
       -        double epsilon, p_f_top, r_norm_max = NAN;
       +        double epsilon, p_f_top, omega, r_norm_max = NAN;
       +        copy_values(sim->p_f_dot, sim->p_f_dot_old, sim->nz);
       +        
        
                /* choose integration method, parameter in [0.0; 1.0]
                 * epsilon = 0.0: explicit
                 * epsilon = 0.5: Crank-Nicolson
                 * epsilon = 1.0: implicit */
                epsilon = 0.5;
       +        omega = 1e-4;
        
                for (i = 0; i < sim->nz; ++i)
                        sim->p_f_dot_impl[i] = sim->p_f_dot_expl[i] = 0.0;
       t@@ -298,7 +301,9 @@ darcy_solver_1d(struct simulation *sim,
                                          + (1.0 - epsilon) * sim->p_f_dot_expl[i];
        
                set_fluid_bcs(sim->p_f_ghost, sim, p_f_top);
       -
       +        for (i = 0; i < sim->nz; ++i)
       +                    sim->p_f_dot[i] = omega * sim->p_f_dot[i]
       +                                      + (1.0 - omega) * sim->p_f_dot_old[i];
        #ifdef DEBUG
                printf(".. epsilon = %g\n", epsilon);
                puts(".. p_f_dot_expl:");
 (DIR) diff --git a/simulation.c b/simulation.c
       t@@ -164,6 +164,7 @@ prepare_arrays(struct simulation *sim)
                sim->old_val = empty(sim->nz);
                sim->fluid_old_val = empty(sim->nz);
                sim->tmp_ghost = empty(sim->nz + 2);
       +        sim->p_f_dot_old = zeros(sim->nz);
        }
        
        void
       t@@ -194,6 +195,7 @@ free_arrays(struct simulation *sim)
                free(sim->old_val);
                free(sim->fluid_old_val);
                free(sim->tmp_ghost);
       +        free(sim->p_f_dot_old);
        }
        
        static void
 (DIR) diff --git a/simulation.h b/simulation.h
       t@@ -125,6 +125,7 @@ struct simulation {
                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);