wip - bddc

examples/petsc/Makefile

@@ -46,6 +46,11 @@ area.o = $(area.c:%.c=$(OBJDIR)/%.o)
 area: $(area.o) | $(PETSc.pc) $(ceed.pc)
 	$(call quiet,LINK.o) $(CEED_LDFLAGS) $^ $(LOADLIBES) $(LDLIBS) -o $@
 
+bddc.c := bddc.c $(utils.c)
+bddc.o = $(bddc.c:%.c=$(OBJDIR)/%.o)
+bddc: $(bddc.o) | $(PETSc.pc) $(ceed.pc)
+	$(call quiet,LINK.o) $(CEED_LDFLAGS) $^ $(LOADLIBES) $(LDLIBS) -o $@
+
 bps.c := bps.c $(utils.c)
 bps.o = $(bps.c:%.c=$(OBJDIR)/%.o)
 bps: $(bps.o) | $(PETSc.pc) $(ceed.pc)
@@ -89,7 +94,7 @@ print: $(PETSc.pc) $(ceed.pc)
 	@true
 
 clean:
-	$(RM) -r $(OBJDIR) *.vtu area bps bpsraw bpssphere multigrid
+	$(RM) -r $(OBJDIR) *.vtu area bddc bps bpsraw bpssphere multigrid
 
 $(PETSc.pc):
 	$(if $(wildcard $@),,$(error \

examples/petsc/bddc.h

@@ -0,0 +1,55 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+// reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+#ifndef bddc_h
+#define bddc_h
+
+#include "include/bpsproblemdata.h"
+#include "include/petscmacros.h"
+#include "include/petscutils.h"
+#include "include/matops.h"
+#include "include/structs.h"
+#include "include/libceedsetup.h"
+
+#include <ceed.h>
+#include <petsc.h>
+#include <petscdmplex.h>
+#include <petscfe.h>
+#include <petscsys.h>
+#include <stdbool.h>
+#include <string.h>
+
+#if PETSC_VERSION_LT(3,12,0)
+#ifdef PETSC_HAVE_CUDA
+#include <petsccuda.h>
+// Note: With PETSc prior to version 3.12.0, providing the source path to
+//       include 'cublas_v2.h' will be needed to use 'petsccuda.h'.
+#endif
+#endif
+
+// -----------------------------------------------------------------------------
+// Command Line Options
+// -----------------------------------------------------------------------------
+
+// Coarsening options
+typedef enum {
+  INJECTION_SCALED = 0, INJECTION_HARMONIC = 1
+} InjectionType;
+static const char *const injection_types [] = {"scaled", "harmonic",
+                                               "InjectionType", "INJECTION", 0
+                                              };
+
+#endif // bddc_h

examples/petsc/include/libceedsetup.h

@@ -17,5 +17,9 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
 PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels,
                                       CeedInt num_comp_u, CeedData *data, CeedInt *leveldegrees,
                                       CeedQFunction qf_restrict, CeedQFunction qf_prolong);
+PetscErrorCode SetupLibceedBDDC(DM dm_vertex, CeedData data_fine,
+                                CeedData data_vertex,
+                                PetscInt g_vertex_size, PetscInt xl_vertex_size,
+                                BPData bp_data);
 
 #endif // libceedsetup_h

examples/petsc/include/matops.h

@@ -13,6 +13,9 @@ PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y);
 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx);
 PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y);
 PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y);
+PetscErrorCode PCShellApply_BDDC(Mat A, Vec X, Vec Y);
+PetscErrorCode PCShellSetup_BDDC(PC pc);
+PetscErrorCode MatMult_BDDCSchur(Mat S, Vec X, Vec Y);
 PetscErrorCode ComputeErrorMax(UserO user, CeedOperator op_error,
                                Vec X, CeedVector target, PetscReal *max_error);
 

examples/petsc/include/petscutils.h

@@ -18,6 +18,8 @@ typedef PetscErrorCode (*BCFunction)(PetscInt dim, PetscReal time,
 PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt num_comp_u,
                                PetscInt topo_dim,
                                bool enforce_bc,  BCFunction bc_func);
+PetscErrorCode SetupVertexDMFromDM(DM dm, DM dm_vertex, PetscInt num_comp_u,
+                                   bool enforce_bc, BCFunction bc_func);
 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P,
     CeedInt topo_dim, CeedInt height, DMLabel domain_label, CeedInt value,
     CeedElemRestriction *elem_restr);

examples/petsc/include/structs.h

@@ -29,6 +29,16 @@ struct UserProlongRestr_ {
   CeedOperator op_prolong, op_restrict;
   Ceed ceed;
 };
+// Data for PETSc PCshell
+typedef struct UserBDDC_ *UserBDDC;
+struct UserBDDC_ {
+  MPI_Comm comm;
+  DM dm, dm_Pi;
+  SNES snes_Pi;
+  KSP ksp_S_Pi;
+  Vec X_loc, Y_loc, X_Pi, Y_Pi, X_Pi_loc, Y_Pi_loc, mult;
+  CeedDataBDDC ceed_data_bddc;
+};
 
 // -----------------------------------------------------------------------------
 // libCEED Data Structs
@@ -45,6 +55,17 @@ struct CeedData_ {
   CeedVector q_data, x_ceed, y_ceed;
 };
 
+// libCEED data struct for BDDC
+typedef struct CeedDataBDDC_ *CeedDataBDDC;
+struct CeedDataBDDC_ {
+  CeedBasis basis_Pi;
+  CeedInt strides[3];
+  CeedElemRestriction elem_restr_Pi, elem_restr_r;
+  CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv, op_inject,
+               op_restrict;
+  CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mult_ceed;
+};
+
 // BP specific data
 typedef struct {
   CeedInt num_comp_x, num_comp_u, topo_dim, q_data_size, q_extra;

examples/petsc/src/libceedsetup.c

@@ -57,11 +57,9 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
   P = degree + 1;
   Q = P + q_extra;
   CeedBasisCreateTensorH1Lagrange(ceed, topo_dim, num_comp_u, P, Q,
-                                  bp_data.q_mode,
-                                  &basis_u);
+                                  bp_data.q_mode, &basis_u);
   CeedBasisCreateTensorH1Lagrange(ceed, topo_dim, num_comp_x, 2, Q,
-                                  bp_data.q_mode,
-                                  &basis_x);
+                                  bp_data.q_mode, &basis_x);
   CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts);
 
   // CEED restrictions
@@ -155,8 +153,7 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
     CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x,
                          CEED_VECTOR_ACTIVE);
     CeedOperatorSetField(op_setup_rhs, "q_data", elem_restr_qd_i,
-                         CEED_BASIS_COLLOCATED,
-                         q_data);
+                         CEED_BASIS_COLLOCATED, q_data);
     CeedOperatorSetField(op_setup_rhs, "true_soln", elem_restr_u_i,
                          CEED_BASIS_COLLOCATED, *target);
     CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u,
@@ -254,4 +251,115 @@ PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels,
   PetscFunctionReturn(0);
 };
 
+// -----------------------------------------------------------------------------
+// Set up libCEED for BDDC interface vertices
+// -----------------------------------------------------------------------------
+PetscErrorCode SetupLibceedBDDC(DM dm_vertex, CeedData data_fine,
+                                CeedDataBDDC data_vertex,
+                                PetscInt g_vertex_size, PetscInt xl_vertex_size,
+BPData bp_data {
+  int ierr;
+  Ceed ceed = data_fine->ceed;
+  CeedBasis basis_Pi, basis_u = data_fine->basis_u;
+  CeedElemRestriction elem_restr_Pi, elem_restr_r;
+  CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv,;
+  CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mask_r_ceed, mask_Gamma_ceed, mask_I_ceed;
+  CeedInt topo_dim, num_comp_u, P, Q, num_qpts, num_elem, elem_size,
+  q_data_size = bp_data.q_data_size;
+
+  // CEED basis
+  CeedBasisGetDimension(basis_u, &topo_dim);
+  CeedBasisGetNumComponents(basis_u, &num_comp_u);
+  CeedBasisGetNumNodes1D(basis_u, &P);
+  elem_size = CeedIntPow(P, topo_dim);
+  CeedBasisGetNumQuadraturePoints1D(basis_u, &Q);
+  CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts);
+  CeedScalar *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d;
+  interp_1d = calloc(2*Q * sizeof(CeedScalar));
+  CeedScalar *temp;
+  CeedBasisGetInterp1D(basis_u, &temp);
+  memcpy(interp_1d, temp, Q * sizeof(CeedScalar));
+  memcpy(&interp_1d[1*Q], temp[(P-1)*Q], Q * sizeof(CeedScalar));
+  grad_1d = calloc(2*Q * sizeof(CeedScalar));
+  CeedBasisGetGrad1D(basis_u, &temp);
+  memcpy(grad_1d, temp, Q * sizeof(CeedScalar));
+  memcpy(&grad_1d[1*Q], temp[(P-1)*Q], Q * sizeof(CeedScalar));
+  q_ref_1d = calloc(Q * sizeof(CeedScalar));
+  CeedBasisGetQRef(basis_u, &temp);
+  memcpy(q_ref_1d, temp, Q * sizeof(CeedScalar));
+  q_weight_1d = calloc(Q * sizeof(CeedScalar));
+  CeedBasisGetQWeights(basis_u, &temp);
+  memcpy(q_weight_1d, temp, Q * sizeof(CeedScalar));
+  CeedBasisCreateTensorH1(ceed, topo_dim, num_comp_u, 1, Q,
+                          interp_1d, grad_1d, q_ref_1d,
+                          q_weight_1d, &basis_Pi);
+
+  // CEED restrictions
+  // -- Interface vertex restriction
+  ierr = CreateRestrictionFromPlex(ceed, dm_vertex, P, topo_dim, 0, 0, 0, &elem_restr_Pi);
+  CHKERRQ(ierr);
+
+  // -- Subdomain restriction
+  ierr = DMPlexGetHeightStratum(dm_vertex, 0, &c_start, &c_end); CHKERRQ(ierr);
+  num_elem = c_end - c_start;
+  CeedInt strides = [num_comp_u, 1, num_comp_u*elem_size];
+  CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, num_comp_u,
+                                   num_comp_u *num_elem*elem_size,
+                                   strides, &elem_restr_r);
+
+  // Create the persistent vectors that will be needed
+  CeedVectorCreate(ceed, xl_vertex_size, &x_Pi_ceed);
+  CeedVectorCreate(ceed, xl_vertex_size, &y_Pi_ceed);
+  CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &x_r_ceed);
+  CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &y_r_ceed);
+
+  // Create the mass or diff operator
+  CeedQFunction qf_apply = data_fine->qf_apply;
+  // -- Interface nodes
+  CeedOperatorSetField(op_Pi_Pi, "u", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_Pi_Pi, "q_data", data_fine->elem_restr_qd_i,
+                       CEED_BASIS_COLLOCATED, data_fine->q_data);
+  CeedOperatorSetField(op_Pi_Pi, "v", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
+  // -- Interface vertices to subdomain
+  CeedOperatorSetField(op_Pi_r, "u", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_Pi_r, "q_data", data_fine->elem_restr_qd_i,
+                       CEED_BASIS_COLLOCATED, data_fine->q_data);
+  CeedOperatorSetField(op_Pi_r, "v", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
+  // -- Subdomain to interface vertices
+  CeedOperatorSetField(op_r_Pi, "u", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_r_Pi, "q_data", data_fine->elem_restr_qd_i,
+                       CEED_BASIS_COLLOCATED, data_fine->q_data);
+  CeedOperatorSetField(op_r_Pi, "v", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
+  // -- Subdomain to subdomain
+  CeedOperatorSetField(op_r_r, "u", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
+  CeedOperatorSetField(op_r_r, "q_data", data_fine->elem_restr_qd_i,
+                       CEED_BASIS_COLLOCATED, data_fine->q_data);
+  CeedOperatorSetField(op_r_r, "v", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
+  // -- Subdomain FDM inverse
+  CeedOperatorCreateFDMElementInverse(op_r_r, &op_r_r_inv, CEED_REQUEST_IMMEDIATE);
+
+  // Injection and restriction operators
+  CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_NONE, CEED_EVAL_NONE,
+                              &qf_inject);
+  CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_NONE, CEED_EVAL_NONE,
+                              &qf_restrict);
+
+  // Save libCEED data required for level
+  data_vertex->basis_Pi = basis_Pi;
+  data_vertex->elem_restr_Pi = elem_restr_Pi;
+  data_vertex->elem_restr_r = elem_restr_r;
+  data_vertex->op_Pi_r = op_Pi_r;
+  data_vertex->op_r_Pi = op_r_Pi;
+  data_vertex->op_Pi_Pi = op_Pi_Pi;
+  data_vertex->op_r_r = op_r_r;
+  data_vertex->op_r_r_inv = op_r_r_inv;
+  data_vertex->x_Pi_ceed = x_Pi_ceed;
+  data_vertex->y_Pi_ceed = y_Pi_ceed;
+  data_vertex->x_r_ceed = x_r_ceed;
+  data_vertex->y_r_ceed = y_r_ceed;
+
+  PetscFunctionReturn(0);
+};
+
+
 // -----------------------------------------------------------------------------