Skip to content

Commit cf2dd0e

Browse files
jeremyltjeremylt
committed
wip
1 parent f5aaafb commit cf2dd0e

File tree

6 files changed

+200
-98
lines changed

6 files changed

+200
-98
lines changed

‎examples/petsc/bddc.c‎

Lines changed: 128 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -46,24 +46,25 @@ int main(int argc, char **argv) {
4646
char filename[PETSC_MAX_PATH_LEN],
4747
ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self";
4848
double my_rt_start, my_rt, rt_min, rt_max;
49-
PetscInt degree = 3, q_extra, l_size, xl_size, g_size, l_vertex_size,
50-
xl_vertex_size, g_vertex_size, dim = 3, m_elem[3] = {3, 3, 3},
51-
num_comp_u = 1;
49+
PetscInt degree = 3, q_extra, l_size, xl_size, g_size, l_Pi_size,
50+
xl_Pi_size, g_Pi_size, dim = 3, m_elem[3] = {3, 3, 3}, num_comp_u = 1;
5251
PetscScalar *r;
5352
PetscScalar eps = 1.0;
5453
PetscBool test_mode, benchmark_mode, read_mesh, write_solution;
5554
PetscLogStage solve_stage;
56-
DM dm, dm_vertex;
55+
DM dm, dm_Pi;
5756
KSP ksp;
5857
PC pc;
59-
Mat mat_O, mat_vertex, mat_pr;
60-
Vec X, X_loc, X_vertex, X_vertex_loc, mult, rhs, rhs_loc;
58+
Mat mat_O, mat_Pi, mat_pr;
59+
Vec X, X_loc, X_Pi, X_Pi_loc, rhs, rhs_loc, mult, mask_r, mask_Gamma,
60+
mask_I;
6161
PetscMemType mem_type;
62-
UserO user_O, user_vertex;
63-
UserProlongRestr user_pr;
62+
UserO user_O, user_Pi;
63+
UserBDDC user_bddc;
6464
Ceed ceed;
65-
CeedData ceed_data, ceed_data_vertex;
66-
CeedVector rhs_ceed, target;
65+
CeedData ceed_data
66+
CeedDataBDDC ceed_data_bddc;
67+
CeedVector rhs_ceed, mult_ceed, target;
6768
CeedQFunction qf_error, qf_restrict, qf_prolong;
6869
CeedOperator op_error;
6970
BPType bp_choice;
@@ -153,7 +154,7 @@ int main(int argc, char **argv) {
153154
dm = dm_dist;
154155
}
155156
}
156-
ierr = DMClone(dm, &dm_vertex); CHKERRQ(ierr);
157+
ierr = DMClone(dm, &dm_Pi); CHKERRQ(ierr);
157158

158159
// Apply Kershaw mesh transformation
159160
ierr = Kershaw(dm, eps); CHKERRQ(ierr);
@@ -179,10 +180,10 @@ int main(int argc, char **argv) {
179180
bp_options[bp_choice].enforce_bc, bp_options[bp_choice].bc_func);
180181
CHKERRQ(ierr);
181182

182-
// Set up vertex DM
183-
ierr = DMClone(dm, &dm_vertex); CHKERRQ(ierr);
184-
ierr = DMSetVecType(dm_vertex, vec_type); CHKERRQ(ierr);
185-
ierr = SetupVertexDMFromDM(dm, dm_vertex, num_comp_u,
183+
// Set up subdomain vertex DM
184+
ierr = DMClone(dm, &dm_Pi); CHKERRQ(ierr);
185+
ierr = DMSetVecType(dm_Pi, vec_type); CHKERRQ(ierr);
186+
ierr = SetupVertexDMFromDM(dm, dm_Pi, num_comp_u,
186187
bp_options[bp_choice].enforce_bc,
187188
bp_options[bp_choice].bc_func);
188189
CHKERRQ(ierr);
@@ -193,11 +194,11 @@ int main(int argc, char **argv) {
193194
ierr = VecGetSize(X, &g_size); CHKERRQ(ierr);
194195
ierr = DMCreateLocalVector(dm, &X_loc); CHKERRQ(ierr);
195196
ierr = VecGetSize(X_loc &xl_size); CHKERRQ(ierr);
196-
ierr = DMCreateGlobalVector(dm_vertex, &X); CHKERRQ(ierr);
197-
ierr = VecGetLocalSize(X_vertex, &l_vertex_size); CHKERRQ(ierr);
198-
ierr = VecGetSize(X_vertex, &g_vertex_size); CHKERRQ(ierr);
199-
ierr = DMCreateLocalVector(dm_vertex, &X_vertex_loc); CHKERRQ(ierr);
200-
ierr = VecGetSize(X_vertex_loc &xl_vertex_size); CHKERRQ(ierr);
197+
ierr = DMCreateGlobalVector(dm_Pi, &X); CHKERRQ(ierr);
198+
ierr = VecGetLocalSize(X_Pi, &l_Pi_size); CHKERRQ(ierr);
199+
ierr = VecGetSize(X_Pi, &g_Pi_size); CHKERRQ(ierr);
200+
ierr = DMCreateLocalVector(dm_Pi, &X_Pi_loc); CHKERRQ(ierr);
201+
ierr = VecGetSize(X_Pi_loc &xl_Pi_size); CHKERRQ(ierr);
201202

202203
// Operator
203204
ierr = PetscMalloc1(1, &user_O); CHKERRQ(ierr);
@@ -210,18 +211,18 @@ int main(int argc, char **argv) {
210211
ierr = MatShellSetVecType(mat_O, vec_type); CHKERRQ(ierr);
211212

212213
// Interface vertex operator
213-
ierr = PetscMalloc1(1, &user_vertex); CHKERRQ(ierr);
214-
ierr = MatCreateShell(comm, l_vertex_size, l_vertex_size, g_vertex_size,
215-
g_vertex_size, user_vertex, &mat_vertex); CHKERRQ(ierr);
214+
ierr = PetscMalloc1(1, &user_Pi); CHKERRQ(ierr);
215+
ierr = MatCreateShell(comm, l_Pi_size, l_Pi_size, g_Pi_size,
216+
g_Pi_size, user_Pi, &mat_Pi); CHKERRQ(ierr);
216217
ierr = MatShellSetOperation(mat_O, MATOP_MULT,
217218
(void(*)(void))MatMult_Ceed); CHKERRQ(ierr);
218-
ierr = MatShellSetOperation(mat_vertex, MATOP_GET_DIAGONAL,
219+
ierr = MatShellSetOperation(mat_Pi, MATOP_GET_DIAGONAL,
219220
(void(*)(void))MatGetDiag); CHKERRQ(ierr);
220-
ierr = MatShellSetVecType(mat_vertex, vec_type); CHKERRQ(ierr);
221+
ierr = MatShellSetVecType(mat_Pi, vec_type); CHKERRQ(ierr);
221222

222223
// Injection operator
223224
ierr = PetscMalloc1(1, &user_pr); CHKERRQ(ierr);
224-
ierr = MatCreateShell(comm, l_size, l_vertex_size, g_size, g_vertex_size,
225+
ierr = MatCreateShell(comm, l_size, l_Pi_size, g_size, g_Pi_size,
225226
user_pr, &mat_pr); CHKERRQ(ierr);
226227
ierr = MatShellSetOperation(mat_pr, MATOP_MULT,
227228
(void(*)(void))MatMult_Inject); CHKERRQ(ierr);
@@ -273,9 +274,9 @@ int main(int argc, char **argv) {
273274
CHKERRQ(ierr);
274275

275276
// Set up libCEED operator on interface vertices
276-
ierr = PetscMalloc1(1, &ceed_data_vertex); CHKERRQ(ierr);
277-
ierr = SetupLibceedBDDC(dm, ceed_data, ceed_data_vertex, g_vertex_size,
278-
xl_vertex_size, bp_options[bp_choice]);
277+
ierr = PetscMalloc1(1, &ceed_data_bddc); CHKERRQ(ierr);
278+
ierr = SetupLibceedBDDC(dm, ceed_data, ceed_data_bddc, g_Pi_size,
279+
xl_Pi_size, bp_options[bp_choice]);
279280
CHKERRQ(ierr);
280281

281282
// Gather RHS
@@ -291,12 +292,6 @@ int main(int argc, char **argv) {
291292
CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_INTERP, CEED_EVAL_NONE,
292293
&qf_prolong);
293294

294-
// Set up libCEED BDDC injection operators
295-
// TODO
296-
ierr = CeedInjectionSetup(ceed, num_comp_u, ceed_data, qf_restrict, qf_prolong);
297-
CHKERRQ(ierr);
298-
// TODO
299-
300295
// Create the error QFunction
301296
CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error,
302297
bp_options[bp_choice].error_loc, &qf_error);
@@ -315,8 +310,96 @@ int main(int argc, char **argv) {
315310
CeedOperatorSetField(op_error, "error", ceed_data[fine_level]->elem_restr_u_i,
316311
CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
317312

318-
// Set up Mat
319-
// User Operator
313+
// PETSc pointwise mult vectors
314+
// -- Calculate multiplicity
315+
{
316+
ierr = VecSet(X_loc, 1.0); CHKERRQ(ierr);
317+
318+
// Local-to-global
319+
ierr = VecZeroEntries(X); CHKERRQ(ierr);
320+
ierr = DMLocalToGlobal(dm, X_loc, ADD_VALUES, X);
321+
CHKERRQ(ierr);
322+
ierr = VecZeroEntries(X_loc); CHKERRQ(ierr);
323+
324+
// Global-to-local
325+
ierr = DMGlobalToLocal(dm, X, INSERT_VALUES, X_loc);
326+
CHKERRQ(ierr);
327+
ierr = VecZeroEntries(X); CHKERRQ(ierr);
328+
329+
// CEED vector
330+
PetscScalar *x;
331+
ierr = VecGetArrayAndMemType(X_loc, &x, &mem_type); CHKERRQ(ierr);
332+
CeedInt len;
333+
CeedVectorGetLength(ceed_data->x_ceed, &len);)
334+
CeedVectorCreate(ceed_data->ceed, len, &mult_ceed);
335+
CeedVectorSetArray(mult_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, x);
336+
337+
// Multiplicity
338+
CeedVector e_vec;
339+
CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, NULL, &e_vec);
340+
CeedVectorSetValue(e_vec, 0.0);
341+
CeedElemRestrictionApply(ceed_data->elem_restr_u, CEED_NOTRANSPOSE, mult_ceed,
342+
e_vec, CEED_REQUEST_IMMEDIATE);
343+
CeedVectorSetValue(mult_ceed, 0.0);
344+
CeedElemRestrictionApply(ceed_data->elem_restr_u, CEED_TRANSPOSE, e_vec,
345+
mult_ceed, CEED_REQUEST_IMMEDIATE);
346+
CeedVectorSyncArray(mult_ceed, MemTypeP2C(mem_type));
347+
CeedVectorDestroy(&e_vec);
348+
349+
// Restore vector
350+
ierr = VecRestoreArrayAndMemType(X_loc, &x); CHKERRQ(ierr);
351+
352+
// Multiplicity scaling
353+
ierr = VecReciprocal(X_loc);
354+
355+
// Copy to Ceed vector
356+
ierr = VecGetArrayAndMemType(X_loc, &x, &mem_type); CHKERRQ(ierr);
357+
CeedVectorSetArray(mult_ceed, MemTypeP2C(mem_type), CEED_COPY_VALUES, x);
358+
ierr = VecRestoreArrayAndMemType(X_loc, &x); CHKERRQ(ierr);
359+
ierr = VecZeroEntries(X_loc); CHKERRQ(ierr);
360+
}
361+
362+
// -- Masks for subdomains
363+
{
364+
CeedInt length_r;
365+
CeedVectorGetLength(ceed_data_bddc->x_r_ceed, &length_r);)
366+
ierr = VecDuplicate(X_loc, &mask_r); CHKERRQ(ierr);
367+
ierr = VecSetSizes(mask_r, length_r, PETSC_DECIDE); CHKERRQ(ierr);
368+
ierr = VecDuplicate(mask_r, &mask_Gamma); CHKERRQ(ierr);
369+
ierr = VecDuplicate(mask_r, &mask_I); CHKERRQ(ierr);
370+
371+
// Set mask contents
372+
CeedScalar *mask_r_array, *mask_Gamma_array, *mask_I_array;
373+
CeedInt num_elem, elem_size;
374+
CeedElemRestrictionGetNumElements(ceed_data_bddc->elem_restr_r, &num_elem);
375+
CeedElemRestrictionGetElementSize(ceed_data_bddc->elem_restr_r, &elem_size);
376+
ierr = VecGetArray(mask_r_ceed, &mask_r_array);
377+
ierr = VecGetArray(mask_Gamma_ceed, &mask_Gamma_array);
378+
ierr = VecGetArray(mask_I_ceed, &mask_I_array);
379+
for (CeedInt e=0; e<num_elem; e++) {
380+
for (CeedInt n=0; n<elem_size; n++) {
381+
PetscBool r = (n % P == 0 || n % P == P-1) &&
382+
((n/P) % P == 0 || (n/P) % P == P-1) &&
383+
((n/(P*P)) % P == 0 || (n/(P*P)) % P == P-1);
384+
PetscBool Gamma = (n % P == 0 || n % P == P-1) ||
385+
((n/P) % P == 0 || (n/P) % P == P-1) ||
386+
((n/(P*P)) % P == 0 || (n/(P*P)) % P == P-1);
387+
PetscBool I = !Gamma;
388+
for (CeedInt c=0; c<num_comp_u; c++) {
389+
CeedInt index = ceed_data_bddc->strides[0]*n + ceed_data_bddc->strides[1]*c +
390+
ceed_data_bddc->strides[2]*e;
391+
mask_r_array[index] = r ? 1.0 : 0.0;
392+
mask_Gamma_array[index] = Gamma ? 1.0 : 0.0;
393+
mask_I_array[index] = I ? 1.0 : 0.0;
394+
}
395+
}
396+
}
397+
ierr = VecRestoreArray(mask_r_ceed, &mask_r_array); CHKERRQ(ierr);
398+
ierr = VecRestoreArray(mask_Gamma_ceed, &mask_Gamma_array); CHKERRQ(ierr);
399+
ierr = VecRestoreArray(mask_I_ceed, &mask_I_array); CHKERRQ(ierr);
400+
}
401+
402+
// Set up MatShell user data
320403
user_O->comm = comm;
321404
user_O->dm = dm;
322405
user_O->X_loc = X_loc;
@@ -326,18 +409,9 @@ int main(int argc, char **argv) {
326409
user_O->op = ceed_data->op_apply;
327410
user_O->ceed = ceed;
328411

329-
// Injection/Restriction Operator
330-
user_pr->comm = comm;
331-
user_pr->dmf = dm;
332-
user_pr->dmc = dm_vertex;
333-
user_pr->loc_vec_c = X_loc;
334-
user_pr->loc_vec_f = user_O->Y_loc;
335-
user_pr->mult_vec = mult;
336-
user_pr->ceed_vec_c = user_O->x_ceed;
337-
user_pr->ceed_vec_f = user_O->y_ceed;
338-
user_pr->op_prolong = ceed_data->op_prolong;
339-
user_pr->op_restrict = ceed_data->op_restrict;
340-
user_pr->ceed = ceed;
412+
// TODO
413+
// Set up PCShell user data
414+
// TODO
341415

342416
// Set up KSP
343417
ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr);
@@ -354,10 +428,10 @@ int main(int argc, char **argv) {
354428
ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
355429
{
356430
ierr = PCSetType(pc, PCSHELL); CHKERRQ(ierr);
431+
ierr = PCShellSetContext(pc, user_bddc); CHKERRQ(ierr);
432+
ierr = PCShellSetApply(pc, (void(*)(void))PCShellApply_BDDC); CHKERRQ(ierr);
357433
// TODO
358-
ierr = PCShellSetContext(pc, ); CHKERRQ(ierr);
359-
ierr = PCShellSetApply(pc, ); CHKERRQ(ierr);
360-
ierr = PCShellSetSetup(pc, ); CHKERRQ(ierr);
434+
//ierr = PCShellSetSetup(pc, ); CHKERRQ(ierr);
361435
// TODO
362436
}
363437

@@ -471,9 +545,9 @@ int main(int argc, char **argv) {
471545
ierr = MatDestroy(&mat_pr); CHKERRQ(ierr);
472546
ierr = PetscFree(user_pr); CHKERRQ(ierr);
473547
ierr = CeedDataDestroy(i, ceed_data); CHKERRQ(ierr);
474-
ierr = CeedDataDestroy(i, ceed_data_vertex); CHKERRQ(ierr);
548+
ierr = CeedDataBDDCDestroy(i, ceed_data_bddc); CHKERRQ(ierr);
475549
ierr = DMDestroy(&dm); CHKERRQ(ierr);
476-
ierr = DMDestroy(&dm_vertex); CHKERRQ(ierr);
550+
ierr = DMDestroy(&dm_Pi); CHKERRQ(ierr);
477551
ierr = VecDestroy(&rhs); CHKERRQ(ierr);
478552
ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr);
479553
ierr = KSPDestroy(&ksp); CHKERRQ(ierr);

‎examples/petsc/include/matops.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -13,8 +13,7 @@ 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 MatMult_Inject(Mat A, Vec X, Vec Y);
17-
PetscErrorCode MatMult_Inject_T(Mat A, Vec X, Vec Y);
16+
PetscErrorCode PCShellApply_BDDC(Mat A, Vec X, Vec Y);
1817
PetscErrorCode ComputeErrorMax(UserO user, CeedOperator op_error,
1918
Vec X, CeedVector target, PetscReal *max_error);
2019

‎examples/petsc/include/structs.h‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -29,15 +29,14 @@ struct UserProlongRestr_ {
2929
CeedOperator op_prolong, op_restrict;
3030
Ceed ceed;
3131
};
32-
33-
// Data for PETSc Inject/Inject_T Matshells
34-
typedef struct UserInject_ *UserInject;
35-
struct UserInject_ {
32+
// Data for PETSc PCshell
33+
typedef struct UserBDDC_ *UserBDDC;
34+
struct UserBDDC_ {
3635
MPI_Comm comm;
37-
DM dmc, dmf;
38-
Vec loc_vec_f, loc_vec_Pi, mult_vec;
39-
CeedVector ceed_vec_f, ceed_vec_Pi, ceed_vec_r;
40-
CeedOperator op_inject, op_inject_t;
36+
DM dm, dm_Pi;
37+
Vec X_loc, Y_loc, diag;
38+
CeedVector x_ceed, y_ceed;
39+
CeedOperator op;
4140
Ceed ceed;
4241
};
4342

@@ -60,9 +59,10 @@ struct CeedData_ {
6059
typedef struct CeedDataBDDC_ *CeedDataBDDC;
6160
struct CeedDataBDDC_ {
6261
CeedBasis basis_Pi;
62+
CeedInt strides[3];
6363
CeedElemRestriction elem_restr_Pi, elem_restr_r;
6464
CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv;
65-
CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mask_r, mask_Gamma, mask_I;
65+
CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mult_ceed;
6666
};
6767

6868
// BP specific data

‎examples/petsc/multigrid.c‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,23 @@ const char help[] = "Solve CEED BPs using p-multigrid with PETSc and DMPlex\n";
4040

4141
#include "bps.h"
4242

43+
// The BDDC example uses vectors in three spaces
44+
//
45+
// Fine mesh: Broken mesh: Vertex mesh:
46+
// x----x----x x----x x----x x x x
47+
// | | | | | | |
48+
// | | | | | | |
49+
// | | | | | | |
50+
// | | | | | | |
51+
// x----x----x x----x x----x x x x
52+
//
53+
// Vectors are organized as follows
54+
// - *_Pi : vector on the vertex mesh
55+
// - *_r : vector on the broken mesh, all points but vertices
56+
// - *_Gamma : vector on the broken mesh, face/vertex/edge points
57+
// - *_I : vector on the broken mesh, interior points
58+
// - * : all other vectors are on the fine mesh
59+
4360
int main(int argc, char **argv) {
4461
PetscInt ierr;
4562
MPI_Comm comm;

0 commit comments

Comments
 (0)