Richards problem converged, and checked with MMS

examples/Hdiv-mixed/conv_test.sh

@@ -32,7 +32,7 @@ declare -A run_flags
     run_flags[pc_type]=svd
     if [[ $dim -eq 2 ]];
     then
-        run_flags[problem]=darcy2d
+        run_flags[problem]=richard2d
         run_flags[dm_plex_dim]=$dim
         run_flags[dm_plex_box_faces]=2,2
     else
@@ -43,8 +43,8 @@ declare -A run_flags
 
 declare -A test_flags
     test_flags[res_start]=2
-    test_flags[res_stride]=1
-    test_flags[res_end]=12
+    test_flags[res_stride]=2
+    test_flags[res_end]=10
 
 file_name=conv_test_result.csv
 

examples/Hdiv-mixed/conv_test_result.csv

@@ -1,12 +1,6 @@
 run,mesh_res,error_u,error_p
-0,2,9.13818,0.09469
-1,3,4.60103,0.05186
-2,4,2.75883,0.03199
-3,5,1.81813,0.02132
-4,6,1.27780,0.01505
-5,7,0.93759,0.01108
-6,8,0.71073,0.00842
-7,9,0.55232,0.00655
-8,10,0.43767,0.00520
-9,11,0.35229,0.00418
-10,12,0.28723,0.00340
+0,2,8.22730,0.00852
+1,4,2.50219,0.00291
+2,6,1.60805,0.00181
+3,8,1.02228,0.00114
+4,10,0.70467,0.00078

examples/Hdiv-mixed/include/setup-solvers.h

@@ -19,10 +19,10 @@ PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data,
                          VecType vec_type, SNES snes, KSP ksp, Vec *U);
 PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U,
                               CeedScalar *l2_error_u, CeedScalar *l2_error_p);
-PetscErrorCode PrintOutput(Ceed ceed,
+PetscErrorCode PrintOutput(Ceed ceed, AppCtx app_ctx, PetscBool has_ts,
                            CeedMemType mem_type_backend,
-                           SNES snes, KSP ksp,
+                           TS ts, SNES snes, KSP ksp,
                            Vec U, CeedScalar l2_error_u,
-                           CeedScalar l2_error_p, AppCtx app_ctx);
+                           CeedScalar l2_error_p);
 
 #endif // setup_solvers_h

examples/Hdiv-mixed/include/setup-ts.h

@@ -20,6 +20,6 @@ PetscErrorCode SetupResidualOperatorCtx_P0(MPI_Comm comm, DM dm, Ceed ceed,
 PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
                                void *ctx_residual_ut);
 PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx,
-                              Vec *U, PetscScalar *f_time, TS *ts);
+                              Vec *U, TS *ts);
 
 #endif // setup_ts_h

examples/Hdiv-mixed/include/structs.h

@@ -18,8 +18,7 @@ struct AppCtx_ {
   // Problem type arguments
   PetscFunctionList problems;
   char              problem_name[PETSC_MAX_PATH_LEN];
-  CeedContextFieldLabel    solution_time_label;
-  CeedScalar        t_final;
+  CeedScalar        t_final, t;
 };
 
 // PETSc operator contexts
@@ -31,8 +30,9 @@ struct OperatorApplyContext_ {
   CeedOperator    op_apply;
   DM              dm;
   Ceed            ceed;
-  CeedScalar      t;
-  CeedContextFieldLabel    solution_time_label, final_time_label;
+  CeedScalar      t, dt;
+  CeedContextFieldLabel solution_time_label, final_time_label,
+                        timestep_label;   ;
 };
 
 // libCEED data struct

examples/Hdiv-mixed/main.c

@@ -136,10 +136,10 @@ int main(int argc, char **argv) {
   //PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
   //                                   bc_pressure) );
 
-
   // ---------------------------------------------------------------------------
   // Setup TSSolve for Richard problem
   // ---------------------------------------------------------------------------
+  TS ts;
   if (problem_data->has_ts) {
     // ---------------------------------------------------------------------------
     // Create global initial conditions
@@ -154,46 +154,49 @@ int main(int argc, char **argv) {
                             ceed_data->ctx_initial_u0,
                             ceed_data->ctx_initial_p0,
                             ceed_data->ctx_residual_ut);
-    VecView(U, PETSC_VIEWER_STDOUT_WORLD);
+    //VecView(U, PETSC_VIEWER_STDOUT_WORLD);
+    // Solve Richards problem
+    PetscCall( VecZeroEntries(ceed_data->ctx_residual_ut->X_loc) );
+    PetscCall( VecZeroEntries(ceed_data->ctx_residual_ut->X_t_loc) );
+    PetscCall( TSSolveRichard(dm, ceed_data, app_ctx,
+                              &U, &ts) );
+    //VecView(U, PETSC_VIEWER_STDOUT_WORLD);
   }
 
+  SNES snes;
+  KSP ksp;
   if (!problem_data->has_ts) {
     // ---------------------------------------------------------------------------
-    // Solve PDE
+    // Setup SNES for Darcy problem
     // ---------------------------------------------------------------------------
     // Create SNES
-    SNES snes;
-    KSP ksp;
     PetscCall( SNESCreate(comm, &snes) );
     PetscCall( SNESGetKSP(snes, &ksp) );
     PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) );
     //VecView(U, PETSC_VIEWER_STDOUT_WORLD);
+  }
 
-    // ---------------------------------------------------------------------------
-    // Compute L2 error of mms problem
-    // ---------------------------------------------------------------------------
-    CeedScalar l2_error_u, l2_error_p;
-    PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, &l2_error_u,
-                              &l2_error_p) );
-
-    // ---------------------------------------------------------------------------
-    // Print output results
-    // ---------------------------------------------------------------------------
-    PetscCall( PrintOutput(ceed, mem_type_backend,
-                           snes, ksp, U, l2_error_u, l2_error_p, app_ctx) );
-
-    // ---------------------------------------------------------------------------
-    // Save solution (paraview)
-    // ---------------------------------------------------------------------------
-    PetscViewer viewer;
+  // ---------------------------------------------------------------------------
+  // Compute L2 error of mms problem
+  // ---------------------------------------------------------------------------
+  CeedScalar l2_error_u, l2_error_p;
+  PetscCall( ComputeL2Error(dm, ceed, ceed_data, U, &l2_error_u,
+                            &l2_error_p) );
 
-    PetscCall( PetscViewerVTKOpen(comm,"solution.vtu",FILE_MODE_WRITE,&viewer) );
-    PetscCall( VecView(U, viewer) );
-    PetscCall( PetscViewerDestroy(&viewer) );
+  // ---------------------------------------------------------------------------
+  // Print output results
+  // ---------------------------------------------------------------------------
+  PetscCall( PrintOutput(ceed, app_ctx, problem_data->has_ts, mem_type_backend,
+                         ts, snes, ksp, U, l2_error_u, l2_error_p) );
 
-    PetscCall( SNESDestroy(&snes) );
+  // ---------------------------------------------------------------------------
+  // Save solution (paraview)
+  // ---------------------------------------------------------------------------
+  PetscViewer viewer;
+  PetscCall( PetscViewerVTKOpen(comm,"solution.vtu",FILE_MODE_WRITE,&viewer) );
+  PetscCall( VecView(U, viewer) );
+  PetscCall( PetscViewerDestroy(&viewer) );
 
-  }
   // ---------------------------------------------------------------------------
   // Free objects
   // ---------------------------------------------------------------------------
@@ -203,6 +206,11 @@ int main(int argc, char **argv) {
   PetscCall( DMDestroy(&dm_u0) );
   PetscCall( DMDestroy(&dm_p0) );
   PetscCall( VecDestroy(&U) );
+  if (problem_data->has_ts) {
+    PetscCall( TSDestroy(&ts) );
+  } else {
+    PetscCall( SNESDestroy(&snes) );
+  }
   PetscCall( CeedDataDestroy(ceed_data, problem_data) );
 
   // -- Function list

examples/Hdiv-mixed/problems/darcy2d.c

@@ -54,7 +54,7 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
   // ------------------------------------------------------
   //              Command line Options
   // ------------------------------------------------------
-  CeedScalar kappa = 1., rho_a0 = 998.2, g = 9.8, alpha_a = 1., b_a = 10.;
+  CeedScalar kappa = 10., rho_a0 = 998.2, g = 9.8, alpha_a = 1., b_a = 10.;
   PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL);
   PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL,
                                 kappa, &kappa, NULL));

examples/Hdiv-mixed/problems/richard2d.c

@@ -21,6 +21,7 @@
 #include "../qfunctions/richard-system2d.h"
 #include "../qfunctions/richard-true2d.h"
 #include "../qfunctions/richard-ics2d.h"
+#include "../qfunctions/darcy-error2d.h"
 //#include "../qfunctions/pressure-boundary2d.h"
 #include "petscsystypes.h"
 
@@ -50,14 +51,12 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
   problem_data->rhs_p0_loc              = RichardRhsP02D_loc;
   problem_data->ics_p                   = RichardICsP2D;
   problem_data->ics_p_loc               = RichardICsP2D_loc;
-  //problem_data->ics_p                   = RichardICsP2D;
-  //problem_data->ics_p_loc               = RichardICsP2D_loc;
-  //problem_data->residual                = RichardSystem2D;
-  //problem_data->residual_loc            = RichardSystem2D_loc;
+  problem_data->residual                = RichardSystem2D;
+  problem_data->residual_loc            = RichardSystem2D_loc;
   //problem_data->jacobian                = JacobianRichardSystem2D;
   //problem_data->jacobian_loc            = JacobianRichardSystem2D_loc;
-  //problem_data->error                   = DarcyError2D;
-  //problem_data->error_loc               = DarcyError2D_loc;
+  problem_data->error                   = DarcyError2D;
+  problem_data->error_loc               = DarcyError2D_loc;
   //problem_data->bc_pressure             = BCPressure2D;
   //problem_data->bc_pressure_loc         = BCPressure2D_loc;
   problem_data->has_ts                  = PETSC_TRUE;
@@ -81,7 +80,7 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
                                 rho_a0, &rho_a0, NULL));
   PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
                                 beta, &beta, NULL));
-  app_ctx->t_final = 5.;
+  app_ctx->t_final = 0.5;
   PetscCall( PetscOptionsScalar("-t_final", "End time", NULL,
                                 app_ctx->t_final, &app_ctx->t_final, NULL));
   PetscOptionsEnd();
@@ -105,10 +104,15 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
                                      offsetof(struct RICHARDContext_, t), 1, "current solver time");
   CeedQFunctionContextRegisterDouble(richard_context, "final_time",
                                      offsetof(struct RICHARDContext_, t_final), 1, "final time");
+  CeedQFunctionContextRegisterDouble(richard_context, "time_step",
+                                     offsetof(struct RICHARDContext_, dt), 1, "time step");
   problem_data->true_qfunction_ctx = richard_context;
   CeedQFunctionContextReferenceCopy(richard_context,
                                     &problem_data->rhs_u0_qfunction_ctx);
-
+  CeedQFunctionContextReferenceCopy(richard_context,
+                                    &problem_data->residual_qfunction_ctx);
+  CeedQFunctionContextReferenceCopy(richard_context,
+                                    &problem_data->error_qfunction_ctx);
   PetscCall( PetscFree(richard_ctx) );
 
   PetscFunctionReturn(0);

examples/Hdiv-mixed/qfunctions/darcy-system2d.h

@@ -172,7 +172,7 @@ CEED_QFUNCTION(JacobianDarcySystem2D)(void *ctx, CeedInt Q,
 
   // *INDENT-ON*
   DARCYContext  context = (DARCYContext)ctx;
-  const CeedScalar  kappa   = context->kappa;
+  const CeedScalar kappa    = context->kappa;
   const CeedScalar rho_a0   = context->rho_a0;
   const CeedScalar g        = context->g;
   const CeedScalar alpha_a  = context->alpha_a;

examples/Hdiv-mixed/qfunctions/darcy-true2d.h

@@ -77,17 +77,20 @@ CEED_QFUNCTION(DarcyTrue2D)(void *ctx, const CeedInt Q,
     CeedScalar x = coords[i+0*Q], y = coords[i+1*Q];  
     CeedScalar psi    = sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
     CeedScalar psi_x  = PI_DOUBLE*cos(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
+    CeedScalar psi_xx  = -PI_DOUBLE*PI_DOUBLE*psi;
     CeedScalar psi_y  = PI_DOUBLE*sin(PI_DOUBLE*x)*cos(PI_DOUBLE*y);
-
+    CeedScalar psi_yy  = -PI_DOUBLE*PI_DOUBLE*psi;
     // k_r = b_a + alpha_a * (1 - x*y)
     CeedScalar k_r = b_a + alpha_a*(1-x*y);
+    CeedScalar k_rx = -alpha_a*y;
+    CeedScalar k_ry = -alpha_a*x;
     // rho = rho_a/rho_a0
     CeedScalar rho = 1.;
     // u = -rho*k_r*K *[grad(\psi) - rho*g_u]
-    CeedScalar u[2] = {-rho*k_r*kappa*psi_x, -rho*k_r*kappa*(psi_y-1)};
-    CeedScalar div_u = -rho*kappa*(-alpha_a*y*psi_x - k_r*PI_DOUBLE*PI_DOUBLE*psi
-                                   -alpha_a*x*psi_y - k_r*PI_DOUBLE*PI_DOUBLE*psi);
-
+    CeedScalar u[2] = {-rho*kappa*k_r*psi_x, 
+                       -rho*kappa*k_r*(psi_y-1)};
+    CeedScalar div_u = -rho*kappa*(k_rx*psi_x + k_r*psi_xx +
+                                   k_ry*(psi_y-1) + k_r*psi_yy);
     // True Force: f = \div(u)
     true_force[i+0*Q] = div_u;
     // True Solution

examples/Hdiv-mixed/qfunctions/darcy-true3d.h

@@ -77,20 +77,26 @@ CEED_QFUNCTION(DarcyTrue3D)(void *ctx, const CeedInt Q,
     CeedScalar x = coords[i+0*Q], y = coords[i+1*Q], z = coords[i+2*Q];
     CeedScalar psi   = sin(PI_DOUBLE*x) * sin(PI_DOUBLE*y) * sin(PI_DOUBLE*z);
     CeedScalar psi_x = PI_DOUBLE*cos(PI_DOUBLE*x) *sin(PI_DOUBLE*y) *sin(PI_DOUBLE*z);
+    CeedScalar psi_xx = -PI_DOUBLE*PI_DOUBLE*psi;
     CeedScalar psi_y = PI_DOUBLE*sin(PI_DOUBLE*x) *cos(PI_DOUBLE*y) *sin(PI_DOUBLE*z);
+    CeedScalar psi_yy = -PI_DOUBLE*PI_DOUBLE*psi;
     CeedScalar psi_z = PI_DOUBLE*sin(PI_DOUBLE*x) *sin(PI_DOUBLE*y) *cos(PI_DOUBLE*z);
+    CeedScalar psi_zz = -PI_DOUBLE*PI_DOUBLE*psi;
 
     // k_r = b_a + alpha_a * (psi - x2)
     CeedScalar k_r = b_a + alpha_a*(1-x*y*z);
+    CeedScalar k_rx = -alpha_a*y*z;
+    CeedScalar k_ry = -alpha_a*x*z;
+    CeedScalar k_rz = -alpha_a*x*y;
     // rho = rho_a/rho_a0
     CeedScalar rho = 1.;
     // u = -rho*k_r*K *[grad(\psi) - rho*g_u]
-    CeedScalar u[3] = {-rho*k_r*kappa*psi_x,
-                       -rho*k_r*kappa*psi_y,
-                       -rho*k_r*kappa*(psi_z-1)};
-    CeedScalar div_u = -rho*kappa*(-alpha_a*y*z*psi_x - k_r*PI_DOUBLE*PI_DOUBLE*psi
-                                   -alpha_a*x*z*psi_y - k_r*PI_DOUBLE*PI_DOUBLE*psi
-                                   -alpha_a*x*y*psi_z - k_r*PI_DOUBLE*PI_DOUBLE*psi);
+    CeedScalar u[3] = {-rho*kappa*k_r*psi_x,
+                       -rho*kappa*k_r*psi_y,
+                       -rho*kappa*k_r*(psi_z-1)};
+    CeedScalar div_u = -rho*kappa*(k_rx*psi_x + k_r*psi_xx +
+                                   k_ry*psi_y + k_r*psi_yy + 
+                                   k_rz*(psi_z-1) + k_r*psi_zz);
 
     // True Force: f = \div(u)
     true_force[i+0*Q] = div_u;

examples/Hdiv-mixed/qfunctions/richard-ics2d.h

@@ -66,7 +66,7 @@ struct RICHARDContext_ {
   CeedScalar rho_a0;
   CeedScalar alpha_a, b_a;
   CeedScalar beta, p0;
-  CeedScalar t, t_final;
+  CeedScalar t, t_final, dt;
   CeedScalar gamma;
 };
 #endif
@@ -105,8 +105,9 @@ CEED_QFUNCTION(RichardRhsU02D)(void *ctx, const CeedInt Q,
     CeedScalar k_r = b_a + alpha_a*(1-x*y);
     // rho = rho_a/rho_a0
     CeedScalar rho = 1.;
-    // ue = -rho*k_r*K *[grad(\psi) - rho*g_u]
-    CeedScalar ue[2] = {-rho*k_r*kappa*psi1_x, -rho*k_r*kappa*(psi1_y-1)};
+    // ue = -rho*k_r*K *[grad(\psi)]
+    CeedScalar ue[2] = {-rho*k_r*kappa*psi1_x,
+                        -rho*k_r*kappa*psi1_y};
     CeedScalar rhs1[2];
     // rhs = (v, ue) = J^T*ue*w
     AlphaMatTransposeVecMult2x2(w[i], J, ue, rhs1);
@@ -184,10 +185,10 @@ CEED_QFUNCTION(RichardRhsP02D)(void *ctx, const CeedInt Q,
                                 {dxdX[0][1][i], dxdX[1][1][i]}};
     const CeedScalar det_J = MatDet2x2(J);
     // psi = exp(-gamma*t)*sin(pi*x)*sin(pi*y)
-    CeedScalar psi_e = sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
+    CeedScalar psi1 = sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
 
     // rhs = (q, pe) = pe*w*det_J
-    rhs_p0[i] = psi_e*w[i]*det_J;
+    rhs_p0[i] = psi1*w[i]*det_J;
   } // End of Quadrature Point Loop
   return 0;
 }