Merge pull request #1510 from CEED/jrwrigh/fluids_matceed

fluids: Add MatCeed

examples/fluids/Makefile

@@ -41,6 +41,7 @@ AFLAGS ?= -fsanitize=address #-fsanitize=undefined -fno-omit-frame-pointer
 CFLAGS += $(if $(ASAN),$(AFLAGS))
 FFLAGS += $(if $(ASAN),$(AFLAGS))
 LDFLAGS += $(if $(ASAN),$(AFLAGS))
+CPPFLAGS += -I./include
 
 OBJDIR := build
 SRCDIR := src

examples/fluids/include/mat-ceed-impl.h

@@ -0,0 +1,69 @@
+#ifndef MAT_CEED_IMPL_H
+#define MAT_CEED_IMPL_H
+
+#include <ceed.h>
+#include <petscdm.h>
+#include <petscmat.h>
+
+#if defined(__clang_analyzer__)
+#define MATCEED_EXTERN extern
+#elif defined(__cplusplus)
+#define MATCEED_EXTERN extern "C"
+#else
+#define MATCEED_EXTERN extern
+#endif
+
+#if defined(__clang_analyzer__)
+#define MATCEED_INTERN
+#else
+#define MATCEED_INTERN MATCEED_EXTERN __attribute__((visibility("hidden")))
+#endif
+
+/**
+  @brief Calls a libCEED function and then checks the resulting error code.
+  If the error code is non-zero, then a PETSc error is set with the libCEED error message.
+**/
+#define PetscCallCeed(ceed_, ...)                                   \
+  do {                                                              \
+    int ierr_q_ = __VA_ARGS__;                                      \
+    if (ierr_q_ != CEED_ERROR_SUCCESS) {                            \
+      const char *error_message;                                    \
+      CeedGetErrorMessage(ceed_, &error_message);                   \
+      SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "%s", error_message); \
+    }                                                               \
+  } while (0)
+
+// MatCeed context for applying composite CeedOperator on a DM
+typedef struct MatCeedContext_private *MatCeedContext;
+struct MatCeedContext_private {
+  Ceed           ceed;
+  char          *name, *internal_mat_type;
+  PetscMemType   mem_type;
+  PetscInt       ref_count, num_mats_assembled_full, num_mats_assembled_pbd;
+  PetscBool      is_destroyed, is_ceed_pbd_valid, is_ceed_vpbd_valid;
+  PetscLogEvent  log_event_mult, log_event_mult_transpose;
+  DM             dm_x, dm_y;
+  Mat           *mats_assembled_full, *mats_assembled_pbd, mat_assembled_full_internal, mat_assembled_pbd_internal;
+  Vec            X_loc, Y_loc_transpose;
+  CeedVector     x_loc, y_loc, coo_values_full, coo_values_pbd;
+  CeedOperator   op_mult, op_mult_transpose;
+  PetscLogDouble flops_mult, flops_mult_transpose;
+};
+
+// Context data
+MATCEED_INTERN PetscErrorCode MatCeedContextCreate(DM dm_x, DM dm_y, Vec X_loc, Vec Y_loc_transpose, CeedOperator op_mult,
+                                                   CeedOperator op_mult_transpose, PetscLogEvent log_event_mult,
+                                                   PetscLogEvent log_event_mult_transpose, MatCeedContext *ctx);
+MATCEED_INTERN PetscErrorCode MatCeedContextReference(MatCeedContext ctx);
+MATCEED_INTERN PetscErrorCode MatCeedContextReferenceCopy(MatCeedContext ctx, MatCeedContext *ctx_copy);
+MATCEED_INTERN PetscErrorCode MatCeedContextDestroy(MatCeedContext ctx);
+
+// Mat Ceed
+MATCEED_INTERN PetscErrorCode MatGetDiagonal_Ceed(Mat A, Vec D);
+MATCEED_INTERN PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y);
+MATCEED_INTERN PetscErrorCode MatMultTranspose_Ceed(Mat A, Vec Y, Vec X);
+
+extern PetscClassId  MATCEED_CLASSID;
+extern PetscLogEvent MATCEED_MULT, MATCEED_MULT_TRANSPOSE;
+
+#endif  // MAT_CEED_IMPL_H

examples/fluids/include/mat-ceed.h

@@ -0,0 +1,41 @@
+#ifndef MAT_CEED_H
+#define MAT_CEED_H
+
+#include <ceed.h>
+#include <petscdm.h>
+#include <petscmat.h>
+
+#define MATCEED "ceed"
+
+#if defined(__clang_analyzer__)
+#define MATCEED_EXTERN extern
+#elif defined(__cplusplus)
+#define MATCEED_EXTERN extern "C"
+#else
+#define MATCEED_EXTERN extern
+#endif
+
+#if defined(__clang_analyzer__)
+#define MATCEED_INTERN
+#else
+#define MATCEED_INTERN MATCEED_EXTERN __attribute__((visibility("hidden")))
+#endif
+
+// Context data
+MATCEED_INTERN PetscErrorCode MatCeedCreate(DM dm_x, DM dm_y, CeedOperator op_mult, CeedOperator op_mult_transpose, Mat *mat);
+MATCEED_INTERN PetscErrorCode MatCeedCopy(Mat mat_ceed, Mat mat_other);
+MATCEED_INTERN PetscErrorCode MatCeedAssembleCOO(Mat mat_ceed, Mat mat_coo);
+MATCEED_INTERN PetscErrorCode MatCeedSetContext(Mat mat, PetscErrorCode (*f)(void *), void *ctx);
+MATCEED_INTERN PetscErrorCode MatCeedGetContext(Mat mat, void *ctx);
+MATCEED_INTERN PetscErrorCode MatCeedSetInnerMatType(Mat mat, MatType type);
+MATCEED_INTERN PetscErrorCode MatCeedGetInnerMatType(Mat mat, MatType *type);
+MATCEED_INTERN PetscErrorCode MatCeedSetOperation(Mat mat, MatOperation op, void (*g)(void));
+MATCEED_INTERN PetscErrorCode MatCeedSetLocalVectors(Mat mat, Vec X_loc, Vec Y_loc_transpose);
+MATCEED_INTERN PetscErrorCode MatCeedGetLocalVectors(Mat mat, Vec *X_loc, Vec *Y_loc_transpose);
+MATCEED_INTERN PetscErrorCode MatCeedRestoreLocalVectors(Mat mat, Vec *X_loc, Vec *Y_loc_transpose);
+MATCEED_INTERN PetscErrorCode MatCeedGetCeedOperators(Mat mat, CeedOperator *op_mult, CeedOperator *op_mult_transpose);
+MATCEED_INTERN PetscErrorCode MatCeedRestoreCeedOperators(Mat mat, CeedOperator *op_mult, CeedOperator *op_mult_transpose);
+MATCEED_INTERN PetscErrorCode MatCeedSetLogEvents(Mat mat, PetscLogEvent log_event_mult, PetscLogEvent log_event_mult_transpose);
+MATCEED_INTERN PetscErrorCode MatCeedGetLogEvents(Mat mat, PetscLogEvent *log_event_mult, PetscLogEvent *log_event_mult_transpose);
+
+#endif  // MAT_CEED_H

examples/fluids/include/petsc_ops.h

@@ -10,6 +10,7 @@
 
 #include <ceed.h>
 #include <petscdm.h>
+#include <petscksp.h>
 
 typedef struct OperatorApplyContext_ *OperatorApplyContext;
 struct OperatorApplyContext_ {
@@ -23,9 +24,6 @@ struct OperatorApplyContext_ {
 PetscErrorCode OperatorApplyContextCreate(DM dm_x, DM dm_y, Ceed ceed, CeedOperator op_apply, CeedVector x_ceed, CeedVector y_ceed, Vec X_loc,
                                           Vec Y_loc, OperatorApplyContext *op_apply_ctx);
 PetscErrorCode OperatorApplyContextDestroy(OperatorApplyContext op_apply_ctx);
-PetscErrorCode MatGetDiag_Ceed(Mat A, Vec D);
-PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y);
-PetscErrorCode CreateMatShell_Ceed(OperatorApplyContext ctx, Mat *mat);
 
 PetscErrorCode DMGetGlobalVectorInfo(DM dm, PetscInt *local_size, PetscInt *global_size, VecType *vec_type);
 PetscErrorCode DMGetLocalVectorInfo(DM dm, PetscInt *local_size, PetscInt *global_size, VecType *vec_type);
@@ -46,4 +44,7 @@ PetscErrorCode ApplyAddCeedOperatorLocalToLocal(Vec X_loc, Vec Y_loc, OperatorAp
 
 PetscErrorCode DMGetLocalVectorInfo(DM dm, PetscInt *local_size, PetscInt *global_size, VecType *vec_type);
 PetscErrorCode DMGetGlobalVectorInfo(DM dm, PetscInt *local_size, PetscInt *global_size, VecType *vec_type);
+
+PetscErrorCode CreateSolveOperatorsFromMatCeed(KSP ksp, Mat mat_ceed, PetscBool assemble, Mat *Amat, Mat *Pmat);
+PetscErrorCode KSPSetFromOptions_WithMatCeed(KSP ksp, Mat mat_ceed);
 #endif  // petsc_ops_h

examples/fluids/navierstokes.c

@@ -23,7 +23,7 @@
 //TESTARGS(name="Blasius, bc_slip") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/blasius.yaml -ts_max_steps 5 -dm_plex_box_faces 3,20,1 -platemesh_nDelta 10 -platemesh_growth 1.2 -bc_outflow 5 -bc_slip 4 -compare_final_state_atol 2E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-bc_slip.bin
 //TESTARGS(name="Blasius, SGS DataDriven Sequential") -ceed {ceed_resource} -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 2e-9 -state_var primitive -ksp_rtol 1e-12 -snes_rtol 1e-12 -stg_mean_only -stg_fluctuating_IC -test_type solver -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin -sgs_model_dd_use_fused false
 //TESTARGS(name="Advection, rotation, cosine") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/advection.yaml -ts_max_steps 0 -advection_ic_type cosine_hill -dm_plex_box_faces 2,1,1 -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-adv-rotation-cosine.bin
-//TESTARGS(name="Gaussian Wave, using MatShell") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 1e-8 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-shell.bin -dm_plex_box_faces 2,2,1 -ts_max_steps 5 -degree 3 -amat_type shell -pmat_pbdiagonal -pc_type vpbjacobi
+//TESTARGS(name="Gaussian Wave, using MatShell") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 1e-8 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-shell.bin -dm_plex_box_faces 2,2,1 -ts_max_steps 5 -degree 3 -amat_type shell -pc_type vpbjacobi
 //TESTARGS(name="Taylor-Green Vortex IC") -ceed {ceed_resource} -problem taylor_green -test_type solver -dm_plex_dim 3 -dm_plex_box_faces 6,6,6 -ts_max_steps 0 -compare_final_state_atol 1e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-taylor-green-IC.bin
 //TESTARGS(name="Blasius, SGS DataDriven Fused") -ceed {ceed_resource} -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 2e-9 -state_var primitive -ksp_rtol 1e-12 -snes_rtol 1e-12 -stg_mean_only -stg_fluctuating_IC -test_type solver -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin
 //TESTARGS(name="Blasius, Anisotropic Differential Filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 5e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_aniso_vandriest.bin -diff_filter_monitor -ts_max_steps 0 -state_var primitive -diff_filter_friction_length 1e-5 -diff_filter_wall_damping_function van_driest -diff_filter_ksp_rtol 1e-8 -diff_filter_grid_based_width -diff_filter_width_scaling 1,0.7,1
@@ -265,8 +265,6 @@ int main(int argc, char **argv) {
   PetscCallCeed(ceed, CeedVectorDestroy(&user->q_ceed));
   PetscCallCeed(ceed, CeedVectorDestroy(&user->q_dot_ceed));
   PetscCallCeed(ceed, CeedVectorDestroy(&user->g_ceed));
-  PetscCallCeed(ceed, CeedVectorDestroy(&user->coo_values_amat));
-  PetscCallCeed(ceed, CeedVectorDestroy(&user->coo_values_pmat));
 
   // -- Bases
   PetscCallCeed(ceed, CeedBasisDestroy(&ceed_data->basis_q));
@@ -316,7 +314,6 @@ int main(int argc, char **argv) {
   PetscCall(OperatorApplyContextDestroy(user->op_rhs_ctx));
   PetscCall(OperatorApplyContextDestroy(user->op_strong_bc_ctx));
   PetscCallCeed(ceed, CeedOperatorDestroy(&user->op_ifunction));
-  PetscCallCeed(ceed, CeedOperatorDestroy(&user->op_ijacobian));
 
   // -- Ceed
   PetscCheck(CeedDestroy(&ceed) == CEED_ERROR_SUCCESS, comm, PETSC_ERR_LIB, "Destroying Ceed object failed");
@@ -340,6 +337,7 @@ int main(int argc, char **argv) {
 
   // -- Matrices
   PetscCall(MatDestroy(&user->interp_viz));
+  PetscCall(MatDestroy(&user->mat_ijacobian));
 
   // -- DM
   PetscCall(DMDestroy(&dm));

examples/fluids/navierstokes.h

@@ -9,6 +9,7 @@
 #define libceed_fluids_examples_navier_stokes_h
 
 #include <ceed.h>
+#include <mat-ceed.h>
 #include <petscts.h>
 #include <stdbool.h>
 
@@ -126,8 +127,7 @@ struct AppCtx_private {
   PetscInt degree;
   PetscInt q_extra;
   // Solver arguments
-  MatType   amat_type;
-  PetscBool pmat_pbdiagonal;
+  MatType amat_type;
   // Post-processing arguments
   PetscInt  checkpoint_interval;
   PetscInt  viz_refine;
@@ -248,11 +248,11 @@ struct User_private {
   Vec                  Q_loc, Q_dot_loc;
   Physics              phys;
   AppCtx               app_ctx;
-  CeedVector           q_ceed, q_dot_ceed, g_ceed, coo_values_amat, coo_values_pmat, x_ceed;
-  CeedOperator         op_rhs_vol, op_ifunction_vol, op_ifunction, op_ijacobian;
+  CeedVector           q_ceed, q_dot_ceed, g_ceed, x_ceed;
+  CeedOperator         op_rhs_vol, op_ifunction_vol, op_ifunction;
+  Mat                  mat_ijacobian;
   KSP                  mass_ksp;
   OperatorApplyContext op_rhs_ctx, op_strong_bc_ctx;
-  bool                 matrices_set_up;
   CeedScalar           time_bc_set;
   SpanStatsData        spanstats;
   NodalProjectionData  grad_velo_proj;

examples/fluids/src/cloptions.c

@@ -102,8 +102,11 @@ PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx, SimpleBC
                                 sizeof(amat_type), &option_set));
     if (option_set) PetscCall(PetscStrallocpy(amat_type, (char **)&app_ctx->amat_type));
   }
-  PetscCall(PetscOptionsBool("-pmat_pbdiagonal", "Assemble only point-block diagonal for Pmat", NULL, app_ctx->pmat_pbdiagonal,
-                             &app_ctx->pmat_pbdiagonal, NULL));
+  {
+    PetscBool option_set;
+    PetscCall(PetscOptionsHasName(NULL, NULL, "-pmat_pbdiagonal", &option_set));
+    if (option_set) PetscCall(PetscPrintf(comm, "Warning! -pmat_pbdiagonal no longer used. Pmat assembly determined from -pc_type setting\n"));
+  }
 
   // Provide default ceed resource if not specified
   if (!ceed_flag) {

examples/fluids/src/differential_filter.c

@@ -78,13 +78,12 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData
   }
 
   {  // Setup LHS Operator and KSP for the differential filtering solve
-    CeedOperator         op_lhs;
-    OperatorApplyContext mat_ctx;
-    Mat                  mat_lhs;
-    CeedInt              num_comp_qd;
-    PetscInt             dim, num_comp_grid_aniso;
-    CeedElemRestriction  elem_restr_grid_aniso;
-    CeedVector           grid_aniso_ceed;
+    CeedOperator        op_lhs;
+    Mat                 mat_lhs;
+    CeedInt             num_comp_qd;
+    PetscInt            dim, num_comp_grid_aniso;
+    CeedElemRestriction elem_restr_grid_aniso;
+    CeedVector          grid_aniso_ceed;
 
     PetscCall(DMGetDimension(user->dm, &dim));
     PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd));
@@ -119,6 +118,7 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData
       }
 
       PetscCallCeed(ceed, CeedQFunctionSetContext(qf_lhs, diff_filter_qfctx));
+      PetscCallCeed(ceed, CeedQFunctionSetUserFlopsEstimate(qf_lhs, 0));
       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "q", num_comp_filter, CEED_EVAL_INTERP));
       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "Grad_q", num_comp_filter * dim, CEED_EVAL_GRAD));
       PetscCallCeed(ceed, CeedQFunctionAddInput(qf_lhs, "anisotropy tensor", num_comp_grid_aniso, CEED_EVAL_NONE));
@@ -149,8 +149,7 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData
       PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_lhs));
       PetscCallCeed(ceed, CeedOperatorDestroy(&op_lhs_sub));
     }
-    PetscCall(OperatorApplyContextCreate(dm_filter, dm_filter, ceed, op_lhs, NULL, NULL, NULL, NULL, &mat_ctx));
-    PetscCall(CreateMatShell_Ceed(mat_ctx, &mat_lhs));
+    PetscCall(MatCeedCreate(dm_filter, dm_filter, op_lhs, NULL, &mat_lhs));
 
     PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm_filter), &diff_filter->ksp));
     PetscCall(KSPSetOptionsPrefix(diff_filter->ksp, "diff_filter_"));
@@ -163,8 +162,7 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData
       PetscCall(KSPSetNormType(diff_filter->ksp, KSP_NORM_NATURAL));
       PetscCall(KSPSetTolerances(diff_filter->ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
     }
-    PetscCall(KSPSetOperators(diff_filter->ksp, mat_lhs, mat_lhs));
-    PetscCall(KSPSetFromOptions(diff_filter->ksp));
+    PetscCall(KSPSetFromOptions_WithMatCeed(diff_filter->ksp, mat_lhs));
 
     PetscCall(MatDestroy(&mat_lhs));
     PetscCallCeed(ceed, CeedOperatorDestroy(&op_lhs));

examples/fluids/src/grid_anisotropy_tensor.c

@@ -73,11 +73,9 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce
   PetscCallCeed(ceed, CeedOperatorSetField(op_mass, "v", *elem_restr_grid_aniso, basis_grid_aniso, CEED_VECTOR_ACTIVE));
 
   {  // -- Setup KSP for L^2 projection
-    Mat                  mat_mass;
-    OperatorApplyContext mass_matop_ctx;
+    Mat mat_mass;
 
-    PetscCall(OperatorApplyContextCreate(grid_aniso_proj->dm, grid_aniso_proj->dm, ceed, op_mass, NULL, NULL, NULL, NULL, &mass_matop_ctx));
-    PetscCall(CreateMatShell_Ceed(mass_matop_ctx, &mat_mass));
+    PetscCall(MatCeedCreate(grid_aniso_proj->dm, grid_aniso_proj->dm, op_mass, NULL, &mat_mass));
 
     PetscCall(KSPCreate(comm, &ksp));
     PetscCall(KSPSetOptionsPrefix(ksp, "grid_anisotropy_tensor_projection_"));
@@ -90,8 +88,7 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce
       PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL));
       PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
     }
-    PetscCall(KSPSetOperators(ksp, mat_mass, mat_mass));
-    PetscCall(KSPSetFromOptions(ksp));
+    PetscCall(KSPSetFromOptions_WithMatCeed(ksp, mat_mass));
     PetscCall(MatDestroy(&mat_mass));
   }