Setup Ifunctiona and TSSolve for Richards problem

Author: rezgarshakeri

examples/Hdiv-mixed/conv_test_result.csv

@@ -1,12 +1,12 @@
 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
+1,3,4.44275,0.05113
+2,4,2.55202,0.03052
+3,5,1.63121,0.01989
+4,6,1.11886,0.01381
+5,7,0.80594,0.01003
+6,8,0.60143,0.00754
+7,9,0.46079,0.00581
+8,10,0.36019,0.00456
+9,11,0.28599,0.00364
+10,12,0.22996,0.00293

examples/Hdiv-mixed/convrate_mixed.png

@@ -6,7 +6,7 @@
 // Convert PETSc MemType to libCEED MemType
 CeedMemType MemTypeP2C(PetscMemType mtype);
 // Destroy libCEED objects
-PetscErrorCode CeedDataDestroy(CeedData ceed_data);
+PetscErrorCode CeedDataDestroy(CeedData ceed_data, ProblemData problem_data);
 // Utility function - essential BC dofs are encoded in closure indices as -(i+1)
 PetscInt Involute(PetscInt i);
 // Utility function to create local CEED restriction from DMPlex
@@ -18,6 +18,6 @@ PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P,
     CeedElemRestriction *elem_restr_oriented, CeedElemRestriction *elem_restr);
 // Set up libCEED for a given degree
 PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
-                            ProblemData problem_data,
-                            PetscInt U_loc_size, CeedData ceed_data);
+                            OperatorApplyContext ctx_residual_ut,
+                            ProblemData problem_data, CeedData ceed_data);
 #endif // setuplibceed_h

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

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

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

@@ -6,12 +6,14 @@
 
 #include "structs.h"
 
-PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *U0);
+PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *Y0,
+                                       OperatorApplyContext ctx_residual_ut);
 PetscErrorCode SetupResidualOperatorCtx_Ut(DM dm, Ceed ceed, CeedData ceed_data,
     OperatorApplyContext ctx_residual_ut);
 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);
+PetscErrorCode TSSolveRichard(DM dm, Ceed ceed, CeedData ceed_data,
+                              AppCtx app_ctx, OperatorApplyContext ctx_residual_ut,
+                              Vec *U, TS *ts);
 
 #endif // setup_ts_h

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

@@ -1,6 +1,8 @@
 #ifndef structs_h
 #define structs_h
 
+#include "ceed/ceed-f64.h"
+#include "petscsystypes.h"
 #include <ceed.h>
 #include <petsc.h>
 
@@ -18,7 +20,7 @@ struct AppCtx_ {
   // Problem type arguments
   PetscFunctionList problems;
   char              problem_name[PETSC_MAX_PATH_LEN];
-  CeedContextFieldLabel    solution_time_label;
+  CeedScalar        t_final;
 };
 
 // 2) richard
@@ -70,9 +72,11 @@ struct OperatorApplyContext_ {
   MPI_Comm        comm;
   Vec             X_loc, Y_loc, X_t_loc;
   CeedVector      x_ceed, y_ceed, x_t_ceed;
-  CeedOperator    op_apply;
+  CeedOperator    op_apply, op_true;
   DM              dm;
   Ceed            ceed;
+  CeedScalar      t;
+  CeedContextFieldLabel    solution_time_label, final_time_label;
 };
 
 // libCEED data struct
@@ -81,10 +85,10 @@ struct CeedData_ {
   CeedBasis            basis_x, basis_u, basis_p, basis_u_face;
   CeedElemRestriction  elem_restr_x, elem_restr_u, elem_restr_U_i,
                        elem_restr_p, elem_restr_p_i;
-  CeedQFunction        qf_residual, qf_jacobian, qf_error, qf_ics;
-  CeedOperator         op_residual, op_jacobian, op_error, op_ics;
-  CeedVector           x_ceed, y_ceed, x_coord, U0_ceed, x_t_ceed;
-  OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut;
+  CeedQFunction        qf_residual, qf_jacobian, qf_error, qf_ics, qf_true;
+  CeedOperator         op_residual, op_jacobian, op_error, op_ics, op_true;
+  CeedVector           x_ceed, y_ceed, x_coord, x_t_ceed;
+  OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error;
 };
 
 // Problem specific data

examples/Hdiv-mixed/include/structs.h

@@ -28,6 +28,8 @@
 //   ./main -pc_type svd -problem darcy3d -dm_plex_filename /path to the mesh file
 //   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
 //   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
+#include "ceed/ceed.h"
+#include <stdio.h>
 const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";
 
 #include "main.h"
@@ -76,6 +78,9 @@ int main(int argc, char **argv) {
   Physics phys_ctx;
   PetscCall( PetscCalloc1(1, &phys_ctx) );
 
+  OperatorApplyContext ctx_residual_ut;
+  PetscCall( PetscCalloc1(1, &ctx_residual_ut) );
+
   // ---------------------------------------------------------------------------
   // Process command line options
   // ---------------------------------------------------------------------------
@@ -113,65 +118,66 @@ int main(int argc, char **argv) {
   // ---------------------------------------------------------------------------
   SetupFE(comm, dm);
 
-  // ---------------------------------------------------------------------------
-  // Create local Force vector
-  // ---------------------------------------------------------------------------
-  Vec U_loc;
-  PetscInt U_loc_size;
-  //CeedVector bc_pressure;
-  PetscCall( DMCreateLocalVector(dm, &U_loc) );
-  // Local size for libCEED
-  PetscCall( VecGetSize(U_loc, &U_loc_size) );
-
   // ---------------------------------------------------------------------------
   // Setup libCEED
   // ---------------------------------------------------------------------------
   // -- Set up libCEED objects
-  PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data,
-                          U_loc_size, ceed_data) );
+  PetscCall( SetupLibceed(dm, ceed, app_ctx, ctx_residual_ut,
+                          problem_data, ceed_data) );
   //CeedVectorView(force_ceed, "%12.8f", stdout);
   //PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
   //                                   bc_pressure) );
 
 
+  // ---------------------------------------------------------------------------
+  // Create Global Solution
+  // ---------------------------------------------------------------------------
+  Vec U; // U = [p,u]
+  PetscCall( DMCreateGlobalVector(dm, &U) );
+
   // ---------------------------------------------------------------------------
   // Setup TSSolve for Richard problem
   // ---------------------------------------------------------------------------
+  TS ts;
+  SNES snes;
+  KSP ksp;
   if (problem_data->has_ts) {
     // ---------------------------------------------------------------------------
     // Create global initial conditions
     // ---------------------------------------------------------------------------
-    Vec U0;
-    CreateInitialConditions(dm, ceed_data, &U0);
-    VecView(U0, PETSC_VIEWER_STDOUT_WORLD);
-    PetscCall( VecDestroy(&U0) );
+    SetupResidualOperatorCtx_Ut(dm, ceed, ceed_data, ctx_residual_ut);
+
+    //SetupResidualOperatorCtx_U0(dm, ceed, ceed_data,
+    //                            ctx_initial);
+    CreateInitialConditions(dm, ceed_data, &U, ctx_residual_ut);
+    //VecView(U, PETSC_VIEWER_STDOUT_WORLD);
+    PetscCall( VecZeroEntries(ctx_residual_ut->X_t_loc) );
+    PetscCall( TSSolveRichard(dm, ceed, ceed_data, app_ctx, ctx_residual_ut,
+                              &U, &ts) );
   }
 
-  // ---------------------------------------------------------------------------
-  // Solve PDE
-  // ---------------------------------------------------------------------------
-  // Create SNES
-  SNES snes;
-  KSP ksp;
-  Vec U;
-  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);
+  if (!problem_data->has_ts) {
+    // ---------------------------------------------------------------------------
+    // Solve PDE
+    // ---------------------------------------------------------------------------
+    // Create SNES
+    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) );
-
+  PetscCall( PrintOutput(ceed, mem_type_backend, ts,
+                         snes, ksp, U, l2_error_u, l2_error_p, app_ctx, problem_data->has_ts) );
   // ---------------------------------------------------------------------------
   // Save solution (paraview)
   // ---------------------------------------------------------------------------
@@ -188,20 +194,23 @@ int main(int argc, char **argv) {
   // Free PETSc objects
   PetscCall( DMDestroy(&dm) );
   PetscCall( VecDestroy(&U) );
-  PetscCall( VecDestroy(&U_loc) );
-  PetscCall( SNESDestroy(&snes) );
+  if (problem_data->has_ts) {
+    PetscCall( TSDestroy(&ts) );
+  } else {
+    PetscCall( SNESDestroy(&snes) );
+  }
 
   // -- Function list
   PetscCall( PetscFunctionListDestroy(&app_ctx->problems) );
 
+  // Free libCEED objects
+  //CeedVectorDestroy(&bc_pressure);
+  PetscCall( CeedDataDestroy(ceed_data, problem_data) );
   // -- Structs
   PetscCall( PetscFree(app_ctx) );
   PetscCall( PetscFree(problem_data) );
   PetscCall( PetscFree(phys_ctx) );
-
-  // Free libCEED objects
-  //CeedVectorDestroy(&bc_pressure);
-  PetscCall( CeedDataDestroy(ceed_data) );
+  PetscCall( PetscFree(ctx_residual_ut) );
   CeedDestroy(&ceed);
 
   return PetscFinalize();

examples/Hdiv-mixed/main.c

@@ -20,6 +20,7 @@
 #include "../include/register-problem.h"
 #include "../qfunctions/richard-system2d.h"
 #include "../qfunctions/richard-true2d.h"
+#include "../qfunctions/darcy-error2d.h"
 #include "../qfunctions/pressure-boundary2d.h"
 #include "petscsystypes.h"
 
@@ -43,12 +44,12 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
   problem_data->ics_loc                 = RichardICs2D_loc;
   problem_data->true_solution           = RichardTrue2D;
   problem_data->true_solution_loc       = RichardTrue2D_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->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->bc_pressure             = BCPressure2D;
   problem_data->bc_pressure_loc         = BCPressure2D_loc;
   problem_data->has_ts                  = PETSC_TRUE;
@@ -72,6 +73,11 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
                                 rho_a0, &rho_a0, NULL));
   PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
                                 beta, &beta, NULL));
+  PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
+                                beta, &beta, NULL));
+  app_ctx->t_final = 1.;
+  PetscCall( PetscOptionsScalar("-t_final", "End time", NULL,
+                                app_ctx->t_final, &app_ctx->t_final, NULL));
   PetscOptionsEnd();
 
   richard_ctx->kappa = kappa;
@@ -82,11 +88,15 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
   richard_ctx->g = g;
   richard_ctx->p0 = p0;
   richard_ctx->gamma = 5.;
+  richard_ctx->t = 0.;
+  richard_ctx->t_final = app_ctx->t_final;
   CeedQFunctionContextCreate(ceed, &richard_context);
   CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES,
                               sizeof(*richard_ctx), richard_ctx);
   CeedQFunctionContextRegisterDouble(richard_context, "time",
                                      offsetof(struct RICHARDContext_, t), 1, "current solver time");
+  CeedQFunctionContextRegisterDouble(richard_context, "final_time",
+                                     offsetof(struct RICHARDContext_, t_final), 1, "final time");
   problem_data->qfunction_context = richard_context;
   CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
                                      FreeContextPetsc);