Get error when try to create ics for pressure

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

@@ -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

examples/Hdiv-mixed/include/structs.h

@@ -38,23 +38,23 @@ struct OperatorApplyContext_ {
 // libCEED data struct
 typedef struct CeedData_ *CeedData;
 struct CeedData_ {
-  CeedBasis            basis_x, basis_u, basis_p, basis_u_face;
+  CeedBasis            basis_x, basis_u, basis_p, basis_u_face, basis_xc;
   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, qf_true;
-  CeedOperator         op_residual, op_jacobian, op_error, op_ics, op_true;
-  CeedVector           x_ceed, y_ceed, x_coord, x0_ceed, x_t_ceed;
+  CeedQFunction        qf_residual, qf_jacobian, qf_error, qf_ics_u, qf_ics_p;
+  CeedOperator         op_residual, op_jacobian, op_error, op_ics_u, op_ics_p;
+  CeedVector           x_ceed, y_ceed, x_coord, u0_ceed, x_t_ceed;
   OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut;
   CeedInt              num_elem;
 };
 
 // Problem specific data
 typedef struct ProblemData_ *ProblemData;
 struct ProblemData_ {
-  CeedQFunctionUser true_solution, residual, jacobian, error, ics,
+  CeedQFunctionUser true_solution, residual, jacobian, error, ics_u, ics_p,
                     bc_pressure;
   const char        *true_solution_loc, *residual_loc, *jacobian_loc,
-        *error_loc, *bc_pressure_loc, *ics_loc;
+        *error_loc, *bc_pressure_loc, *ics_u_loc, *ics_p_loc;
   CeedQuadMode      quadrature_mode;
   CeedInt           elem_node, dim, q_data_size_face;
   CeedQFunctionContext qfunction_context;

examples/Hdiv-mixed/main.c

@@ -143,8 +143,8 @@ int main(int argc, char **argv) {
     // ---------------------------------------------------------------------------
     // Create global initial conditions
     // ---------------------------------------------------------------------------
-    CreateInitialConditions(dm, ceed_data, U_loc, U);
-    VecView(U, PETSC_VIEWER_STDOUT_WORLD);
+    //CreateInitialConditions(dm, ceed_data, U_loc, U);
+    //VecView(U, PETSC_VIEWER_STDOUT_WORLD);
   }
 
   if (!problem_data->has_ts) {
@@ -192,6 +192,7 @@ int main(int argc, char **argv) {
   PetscCall( DMDestroy(&dm) );
   PetscCall( VecDestroy(&U) );
   PetscCall( VecDestroy(&U_loc) );
+  PetscCall( CeedDataDestroy(ceed_data, problem_data) );
 
   // -- Function list
   PetscCall( PetscFunctionListDestroy(&app_ctx->problems) );
@@ -203,9 +204,6 @@ int main(int argc, char **argv) {
 
   // Free libCEED objects
   //CeedVectorDestroy(&bc_pressure);
-  //PetscCall( CeedDataDestroy(ceed_data) );
-  CeedQFunctionDestroy(&ceed_data->qf_true);
-  CeedOperatorDestroy(&ceed_data->op_true);
   CeedDestroy(&ceed);
 
   return PetscFinalize();

examples/Hdiv-mixed/problems/richard2d.c

@@ -39,8 +39,10 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
   problem_data->elem_node               = 4;
   problem_data->q_data_size_face        = 3;
   problem_data->quadrature_mode         = CEED_GAUSS;
-  problem_data->ics                     = RichardICs2D;
-  problem_data->ics_loc                 = RichardICs2D_loc;
+  problem_data->ics_u                   = RichardICsU2D;
+  problem_data->ics_u_loc               = RichardICsU2D_loc;
+  problem_data->ics_p                   = RichardICsP2D;
+  problem_data->ics_p_loc               = RichardICsP2D_loc;
   problem_data->true_solution           = RichardTrue2D;
   problem_data->true_solution_loc       = RichardTrue2D_loc;
   //problem_data->residual                = RichardSystem2D;
@@ -97,6 +99,7 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
   CeedQFunctionContextRegisterDouble(richard_context, "final_time",
                                      offsetof(struct RICHARDContext_, t_final), 1, "final time");
   problem_data->qfunction_context = richard_context;
+  PetscCall( PetscFree(richard_ctx) );
 
   PetscFunctionReturn(0);
 }

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

@@ -21,6 +21,7 @@
 #define DARCY_ERROR2D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 // -----------------------------------------------------------------------------

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

@@ -21,6 +21,7 @@
 #define DARCY_ERROR3D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 

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

@@ -21,7 +21,7 @@
 #define DARCY_SYSTEM2D_H
 
 #include <math.h>
-#include "ceed/ceed-f64.h"
+#include <ceed.h>
 #include "utils.h"
 
 // -----------------------------------------------------------------------------

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

@@ -21,6 +21,7 @@
 #define DARCY_SYSTEM3D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 // -----------------------------------------------------------------------------

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

@@ -21,6 +21,7 @@
 #define DARCY_TRUE2D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 // -----------------------------------------------------------------------------

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

@@ -21,6 +21,7 @@
 #define DARCY_TRUE3D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 // -----------------------------------------------------------------------------

examples/Hdiv-mixed/qfunctions/pressure-boundary2d.h

@@ -11,6 +11,9 @@
 #ifndef pressure_bc_2d_h
 #define pressure_bc_2d_h
 
+#include <math.h>
+#include <ceed.h>
+#include "utils.h"
 // -----------------------------------------------------------------------------
 // Strong form:
 //  u       = -\grad(p)      on \Omega

examples/Hdiv-mixed/qfunctions/pressure-boundary3d.h

@@ -11,6 +11,9 @@
 #ifndef pressure_bc_3d_h
 #define pressure_bc_3d_h
 
+#include <math.h>
+#include <ceed.h>
+#include "utils.h"
 // -----------------------------------------------------------------------------
 // Strong form:
 //  u       = -\grad(p)      on \Omega

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

@@ -21,6 +21,7 @@
 #define RICHARD_SYSTEM2D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 // See Matthew Farthing, Christopher Kees, Cass Miller (2003)

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

@@ -21,6 +21,7 @@
 #define RICHARD_TRUE2D_H
 
 #include <math.h>
+#include <ceed.h>
 #include "utils.h"
 
 // See Matthew Farthing, Christopher Kees, Cass Miller (2003)
@@ -80,15 +81,15 @@ struct RICHARDContext_ {
 // -----------------------------------------------------------------------------
 // Initial conditions for Richard problem
 // -----------------------------------------------------------------------------
-CEED_QFUNCTION(RichardICs2D)(void *ctx, const CeedInt Q,
-                             const CeedScalar *const *in,
-                             CeedScalar *const *out) {
+CEED_QFUNCTION(RichardICsU2D)(void *ctx, const CeedInt Q,
+                              const CeedScalar *const *in,
+                              CeedScalar *const *out) {
   // *INDENT-OFF*
   // Inputs
   const CeedScalar (*coords) = in[0],
                    (*dxdX)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[1];
   // Outputs
-  CeedScalar (*u_0) = out[0], (*psi_0) = out[1];
+  CeedScalar (*u_0) = out[0];
   RICHARDContext  context = (RICHARDContext)ctx;
   const CeedScalar kappa    = context->kappa;
   const CeedScalar alpha_a  = context->alpha_a;
@@ -108,9 +109,6 @@ CEED_QFUNCTION(RichardICs2D)(void *ctx, const CeedInt Q,
   CeedScalar x1 = coords[1+0*Q], y1 = coords[1+1*Q];
   CeedScalar x2 = coords[2+0*Q], y2 = coords[2+1*Q];
   CeedScalar x3 = coords[3+0*Q], y3 = coords[3+1*Q];
-  CeedScalar xc = (x0+x1)/2., yc = (y0+y2)/2.;
-
-  CeedScalar psic = sin(PI_DOUBLE*xc)*sin(PI_DOUBLE*yc);
 
   CeedScalar psi0_x = PI_DOUBLE*cos(PI_DOUBLE*x0)*sin(PI_DOUBLE*y0);
   CeedScalar psi0_y = PI_DOUBLE*sin(PI_DOUBLE*x0)*cos(PI_DOUBLE*y0);
@@ -149,6 +147,7 @@ CEED_QFUNCTION(RichardICs2D)(void *ctx, const CeedInt Q,
   d5 = ue3[0]*nt3[0]+ue3[1]*nt3[1];
   d6 = ue0[0]*nl0[0]+ue0[1]*nl0[1];
   d7 = ue2[0]*nl2[0]+ue2[1]*nl2[1];
+  // 2nd output 
   u_0[0] = d0;
   u_0[1] = d1;
   u_0[2] = d2;
@@ -157,11 +156,31 @@ CEED_QFUNCTION(RichardICs2D)(void *ctx, const CeedInt Q,
   u_0[5] = d5;
   u_0[6] = d6;
   u_0[7] = d7;
-  // 2nd eq
-  psi_0[0] = psic;
-  psi_0[1] = psic;
-  psi_0[2] = psic;
-  psi_0[3] = psic;
+  } // End of Quadrature Point Loop
+  return 0;
+}
+
+CEED_QFUNCTION(RichardICsP2D)(void *ctx, const CeedInt Q,
+                             const CeedScalar *const *in,
+                             CeedScalar *const *out) {
+  // *INDENT-OFF*
+  // Inputs
+  const CeedScalar (*coords) = in[0];
+  // Outputs
+  CeedScalar (*psi_0) = out[0];
+  // Quadrature Point Loop
+  {
+  CeedScalar x0 = coords[0+0*Q], y0 = coords[0+1*Q];
+  CeedScalar x1 = coords[1+0*Q];
+  CeedScalar y2 = coords[2+1*Q];
+  CeedScalar xc = (x0+x1)/2., yc = (y0+y2)/2.;
+  CeedScalar psi_c = sin(PI_DOUBLE*xc)*sin(PI_DOUBLE*yc);
+
+  // 1st output
+  psi_0[0] = psi_c;
+  //psi_0[1] = psi_c;
+  //psi_0[2] = psi_c;
+  //psi_0[3] = psi_c;
   } // End of Quadrature Point Loop
   return 0;
 }