Skip to content

Commit c2b2f48

Browse files
jeremyltjeremylt
committed
wip - bddc
1 parent d8d9c82 commit c2b2f48

File tree

10 files changed

+1047
-7
lines changed

10 files changed

+1047
-7
lines changed

‎examples/petsc/Makefile‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,11 @@ area.o = $(area.c:%.c=$(OBJDIR)/%.o)
4646
area: $(area.o) | $(PETSc.pc) $(ceed.pc)
4747
$(call quiet,LINK.o) $(CEED_LDFLAGS) $^ $(LOADLIBES) $(LDLIBS) -o $@
4848

49+
bddc.c := bddc.c $(utils.c)
50+
bddc.o = $(bddc.c:%.c=$(OBJDIR)/%.o)
51+
bddc: $(bddc.o) | $(PETSc.pc) $(ceed.pc)
52+
$(call quiet,LINK.o) $(CEED_LDFLAGS) $^ $(LOADLIBES) $(LDLIBS) -o $@
53+
4954
bps.c := bps.c $(utils.c)
5055
bps.o = $(bps.c:%.c=$(OBJDIR)/%.o)
5156
bps: $(bps.o) | $(PETSc.pc) $(ceed.pc)
@@ -89,7 +94,7 @@ print: $(PETSc.pc) $(ceed.pc)
8994
@true
9095

9196
clean:
92-
$(RM) -r $(OBJDIR) *.vtu area bps bpsraw bpssphere multigrid
97+
$(RM) -r $(OBJDIR) *.vtu area bddc bps bpsraw bpssphere multigrid
9398

9499
$(PETSc.pc):
95100
$(if $(wildcard $@),,$(error \

‎examples/petsc/bddc.c‎

Lines changed: 610 additions & 0 deletions
Large diffs are not rendered by default.

‎examples/petsc/bddc.h‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2+
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3+
// reserved. See files LICENSE and NOTICE for details.
4+
//
5+
// This file is part of CEED, a collection of benchmarks, miniapps, software
6+
// libraries and APIs for efficient high-order finite element and spectral
7+
// element discretizations for exascale applications. For more information and
8+
// source code availability see http://github.com/ceed.
9+
//
10+
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11+
// a collaborative effort of two U.S. Department of Energy organizations (Office
12+
// of Science and the National Nuclear Security Administration) responsible for
13+
// the planning and preparation of a capable exascale ecosystem, including
14+
// software, applications, hardware, advanced system engineering and early
15+
// testbed platforms, in support of the nation's exascale computing imperative.
16+
17+
#ifndef bddc_h
18+
#define bddc_h
19+
20+
#include "include/bpsproblemdata.h"
21+
#include "include/petscmacros.h"
22+
#include "include/petscutils.h"
23+
#include "include/matops.h"
24+
#include "include/structs.h"
25+
#include "include/libceedsetup.h"
26+
27+
#include <ceed.h>
28+
#include <petsc.h>
29+
#include <petscdmplex.h>
30+
#include <petscfe.h>
31+
#include <petscsys.h>
32+
#include <stdbool.h>
33+
#include <string.h>
34+
35+
#if PETSC_VERSION_LT(3,12,0)
36+
#ifdef PETSC_HAVE_CUDA
37+
#include <petsccuda.h>
38+
// Note: With PETSc prior to version 3.12.0, providing the source path to
39+
// include 'cublas_v2.h' will be needed to use 'petsccuda.h'.
40+
#endif
41+
#endif
42+
43+
// -----------------------------------------------------------------------------
44+
// Command Line Options
45+
// -----------------------------------------------------------------------------
46+
47+
// Coarsening options
48+
typedef enum {
49+
INJECTION_SCALED = 0, INJECTION_HARMONIC = 1
50+
} InjectionType;
51+
static const char *const injection_types [] = {"scaled", "harmonic",
52+
"InjectionType", "INJECTION", 0
53+
};
54+
55+
#endif // bddc_h

‎examples/petsc/include/libceedsetup.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,5 +17,9 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
1717
PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels,
1818
CeedInt num_comp_u, CeedData *data, CeedInt *leveldegrees,
1919
CeedQFunction qf_restrict, CeedQFunction qf_prolong);
20+
PetscErrorCode SetupLibceedBDDC(DM dm_vertex, CeedData data_fine,
21+
CeedData data_vertex,
22+
PetscInt g_vertex_size, PetscInt xl_vertex_size,
23+
BPData bp_data);
2024

2125
#endif // libceedsetup_h

‎examples/petsc/include/matops.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,9 @@ PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y);
1313
PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx);
1414
PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y);
1515
PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y);
16+
PetscErrorCode PCShellApply_BDDC(Mat A, Vec X, Vec Y);
17+
PetscErrorCode PCShellSetup_BDDC(PC pc);
18+
PetscErrorCode MatMult_BDDCSchur(Mat S, Vec X, Vec Y);
1619
PetscErrorCode ComputeErrorMax(UserO user, CeedOperator op_error,
1720
Vec X, CeedVector target, PetscReal *max_error);
1821

‎examples/petsc/include/petscutils.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,8 @@ typedef PetscErrorCode (*BCFunction)(PetscInt dim, PetscReal time,
1818
PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt num_comp_u,
1919
PetscInt topo_dim,
2020
bool enforce_bc, BCFunction bc_func);
21+
PetscErrorCode SetupVertexDMFromDM(DM dm, DM dm_vertex, PetscInt num_comp_u,
22+
bool enforce_bc, BCFunction bc_func);
2123
PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P,
2224
CeedInt topo_dim, CeedInt height, DMLabel domain_label, CeedInt value,
2325
CeedElemRestriction *elem_restr);

‎examples/petsc/include/structs.h‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,16 @@ struct UserProlongRestr_ {
2929
CeedOperator op_prolong, op_restrict;
3030
Ceed ceed;
3131
};
32+
// Data for PETSc PCshell
33+
typedef struct UserBDDC_ *UserBDDC;
34+
struct UserBDDC_ {
35+
MPI_Comm comm;
36+
DM dm, dm_Pi;
37+
SNES snes_Pi;
38+
KSP ksp_S_Pi;
39+
Vec X_loc, Y_loc, X_Pi, Y_Pi, X_Pi_loc, Y_Pi_loc, mult;
40+
CeedDataBDDC ceed_data_bddc;
41+
};
3242

3343
// -----------------------------------------------------------------------------
3444
// libCEED Data Structs
@@ -45,6 +55,17 @@ struct CeedData_ {
4555
CeedVector q_data, x_ceed, y_ceed;
4656
};
4757

58+
// libCEED data struct for BDDC
59+
typedef struct CeedDataBDDC_ *CeedDataBDDC;
60+
struct CeedDataBDDC_ {
61+
CeedBasis basis_Pi;
62+
CeedInt strides[3];
63+
CeedElemRestriction elem_restr_Pi, elem_restr_r;
64+
CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv, op_inject,
65+
op_restrict;
66+
CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mult_ceed;
67+
};
68+
4869
// BP specific data
4970
typedef struct {
5071
CeedInt num_comp_x, num_comp_u, topo_dim, q_data_size, q_extra;

‎examples/petsc/src/libceedsetup.c‎

Lines changed: 114 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -57,11 +57,9 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
5757
P = degree + 1;
5858
Q = P + q_extra;
5959
CeedBasisCreateTensorH1Lagrange(ceed, topo_dim, num_comp_u, P, Q,
60-
bp_data.q_mode,
61-
&basis_u);
60+
bp_data.q_mode, &basis_u);
6261
CeedBasisCreateTensorH1Lagrange(ceed, topo_dim, num_comp_x, 2, Q,
63-
bp_data.q_mode,
64-
&basis_x);
62+
bp_data.q_mode, &basis_x);
6563
CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts);
6664

6765
// CEED restrictions
@@ -155,8 +153,7 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
155153
CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x,
156154
CEED_VECTOR_ACTIVE);
157155
CeedOperatorSetField(op_setup_rhs, "q_data", elem_restr_qd_i,
158-
CEED_BASIS_COLLOCATED,
159-
q_data);
156+
CEED_BASIS_COLLOCATED, q_data);
160157
CeedOperatorSetField(op_setup_rhs, "true_soln", elem_restr_u_i,
161158
CEED_BASIS_COLLOCATED, *target);
162159
CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u,
@@ -254,4 +251,115 @@ PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels,
254251
PetscFunctionReturn(0);
255252
};
256253

254+
// -----------------------------------------------------------------------------
255+
// Set up libCEED for BDDC interface vertices
256+
// -----------------------------------------------------------------------------
257+
PetscErrorCode SetupLibceedBDDC(DM dm_vertex, CeedData data_fine,
258+
CeedDataBDDC data_vertex,
259+
PetscInt g_vertex_size, PetscInt xl_vertex_size,
260+
BPData bp_data {
261+
int ierr;
262+
Ceed ceed = data_fine->ceed;
263+
CeedBasis basis_Pi, basis_u = data_fine->basis_u;
264+
CeedElemRestriction elem_restr_Pi, elem_restr_r;
265+
CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv,;
266+
CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mask_r_ceed, mask_Gamma_ceed, mask_I_ceed;
267+
CeedInt topo_dim, num_comp_u, P, Q, num_qpts, num_elem, elem_size,
268+
q_data_size = bp_data.q_data_size;
269+
270+
// CEED basis
271+
CeedBasisGetDimension(basis_u, &topo_dim);
272+
CeedBasisGetNumComponents(basis_u, &num_comp_u);
273+
CeedBasisGetNumNodes1D(basis_u, &P);
274+
elem_size = CeedIntPow(P, topo_dim);
275+
CeedBasisGetNumQuadraturePoints1D(basis_u, &Q);
276+
CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts);
277+
CeedScalar *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d;
278+
interp_1d = calloc(2*Q * sizeof(CeedScalar));
279+
CeedScalar *temp;
280+
CeedBasisGetInterp1D(basis_u, &temp);
281+
memcpy(interp_1d, temp, Q * sizeof(CeedScalar));
282+
memcpy(&interp_1d[1*Q], temp[(P-1)*Q], Q * sizeof(CeedScalar));
283+
grad_1d = calloc(2*Q * sizeof(CeedScalar));
284+
CeedBasisGetGrad1D(basis_u, &temp);
285+
memcpy(grad_1d, temp, Q * sizeof(CeedScalar));
286+
memcpy(&grad_1d[1*Q], temp[(P-1)*Q], Q * sizeof(CeedScalar));
287+
q_ref_1d = calloc(Q * sizeof(CeedScalar));
288+
CeedBasisGetQRef(basis_u, &temp);
289+
memcpy(q_ref_1d, temp, Q * sizeof(CeedScalar));
290+
q_weight_1d = calloc(Q * sizeof(CeedScalar));
291+
CeedBasisGetQWeights(basis_u, &temp);
292+
memcpy(q_weight_1d, temp, Q * sizeof(CeedScalar));
293+
CeedBasisCreateTensorH1(ceed, topo_dim, num_comp_u, 1, Q,
294+
interp_1d, grad_1d, q_ref_1d,
295+
q_weight_1d, &basis_Pi);
296+
297+
// CEED restrictions
298+
// -- Interface vertex restriction
299+
ierr = CreateRestrictionFromPlex(ceed, dm_vertex, P, topo_dim, 0, 0, 0, &elem_restr_Pi);
300+
CHKERRQ(ierr);
301+
302+
// -- Subdomain restriction
303+
ierr = DMPlexGetHeightStratum(dm_vertex, 0, &c_start, &c_end); CHKERRQ(ierr);
304+
num_elem = c_end - c_start;
305+
CeedInt strides = [num_comp_u, 1, num_comp_u*elem_size];
306+
CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, num_comp_u,
307+
num_comp_u *num_elem*elem_size,
308+
strides, &elem_restr_r);
309+
310+
// Create the persistent vectors that will be needed
311+
CeedVectorCreate(ceed, xl_vertex_size, &x_Pi_ceed);
312+
CeedVectorCreate(ceed, xl_vertex_size, &y_Pi_ceed);
313+
CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &x_r_ceed);
314+
CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &y_r_ceed);
315+
316+
// Create the mass or diff operator
317+
CeedQFunction qf_apply = data_fine->qf_apply;
318+
// -- Interface nodes
319+
CeedOperatorSetField(op_Pi_Pi, "u", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
320+
CeedOperatorSetField(op_Pi_Pi, "q_data", data_fine->elem_restr_qd_i,
321+
CEED_BASIS_COLLOCATED, data_fine->q_data);
322+
CeedOperatorSetField(op_Pi_Pi, "v", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
323+
// -- Interface vertices to subdomain
324+
CeedOperatorSetField(op_Pi_r, "u", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
325+
CeedOperatorSetField(op_Pi_r, "q_data", data_fine->elem_restr_qd_i,
326+
CEED_BASIS_COLLOCATED, data_fine->q_data);
327+
CeedOperatorSetField(op_Pi_r, "v", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
328+
// -- Subdomain to interface vertices
329+
CeedOperatorSetField(op_r_Pi, "u", elem_restr_Pi, basis_u, CEED_VECTOR_ACTIVE);
330+
CeedOperatorSetField(op_r_Pi, "q_data", data_fine->elem_restr_qd_i,
331+
CEED_BASIS_COLLOCATED, data_fine->q_data);
332+
CeedOperatorSetField(op_r_Pi, "v", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
333+
// -- Subdomain to subdomain
334+
CeedOperatorSetField(op_r_r, "u", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
335+
CeedOperatorSetField(op_r_r, "q_data", data_fine->elem_restr_qd_i,
336+
CEED_BASIS_COLLOCATED, data_fine->q_data);
337+
CeedOperatorSetField(op_r_r, "v", elem_restr_r, basis_u, CEED_VECTOR_ACTIVE);
338+
// -- Subdomain FDM inverse
339+
CeedOperatorCreateFDMElementInverse(op_r_r, &op_r_r_inv, CEED_REQUEST_IMMEDIATE);
340+
341+
// Injection and restriction operators
342+
CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_NONE, CEED_EVAL_NONE,
343+
&qf_inject);
344+
CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_NONE, CEED_EVAL_NONE,
345+
&qf_restrict);
346+
347+
// Save libCEED data required for level
348+
data_vertex->basis_Pi = basis_Pi;
349+
data_vertex->elem_restr_Pi = elem_restr_Pi;
350+
data_vertex->elem_restr_r = elem_restr_r;
351+
data_vertex->op_Pi_r = op_Pi_r;
352+
data_vertex->op_r_Pi = op_r_Pi;
353+
data_vertex->op_Pi_Pi = op_Pi_Pi;
354+
data_vertex->op_r_r = op_r_r;
355+
data_vertex->op_r_r_inv = op_r_r_inv;
356+
data_vertex->x_Pi_ceed = x_Pi_ceed;
357+
data_vertex->y_Pi_ceed = y_Pi_ceed;
358+
data_vertex->x_r_ceed = x_r_ceed;
359+
data_vertex->y_r_ceed = y_r_ceed;
360+
361+
PetscFunctionReturn(0);
362+
};
363+
364+
257365
// -----------------------------------------------------------------------------

0 commit comments

Comments
 (0)