@@ -29,6 +29,35 @@ PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) {
2929 PetscFunctionReturn (0 );
3030};
3131
32+ // -----------------------------------------------------------------------------
33+ // Destroy libCEED BDDC objects
34+ // -----------------------------------------------------------------------------
35+ PetscErrorCode CeedDataBDDCDestroy (CeedDataBDDC data ) {
36+ int ierr ;
37+
38+ CeedBasisDestroy (& data -> basis_Pi );
39+ CeedElemRestrictionDestroy (& data -> elem_restr_Pi );
40+ CeedElemRestrictionDestroy (& data -> elem_restr_Pi_r );
41+ CeedElemRestrictionDestroy (& data -> elem_restr_r );
42+ CeedOperatorDestroy (& data -> op_Pi_r ;);
43+ CeedOperatorDestroy (& data -> op_r_Pi );
44+ CeedOperatorDestroy (& data -> op_Pi_Pi );
45+ CeedOperatorDestroy (& data -> op_r_r );
46+ CeedOperatorDestroy (& data -> op_r_r_inv );
47+ CeedOperatorDestroy (& data -> op_inject_Pi );
48+ CeedOperatorDestroy (& data -> op_inject_r );
49+ CeedOperatorDestroy (& data -> op_restrict_Pi );
50+ CeedOperatorDestroy (& data -> op_restrict_r );
51+ CeedVectorDestroy (& data -> x_Pi_ceed );
52+ CeedVectorDestroy (& data -> y_Pi_ceed );
53+ CeedVectorDestroy (& data -> x_r_ceed );
54+ CeedVectorDestroy (& data -> y_r_ceed );
55+ CeedVectorDestroy (& data -> mult_ceed );
56+ ierr = PetscFree (data ); CHKERRQ (ierr );
57+
58+ PetscFunctionReturn (0 );
59+ };
60+
3261// -----------------------------------------------------------------------------
3362// Set up libCEED for a given degree
3463// -----------------------------------------------------------------------------
@@ -57,11 +86,9 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
5786 P = degree + 1 ;
5887 Q = P + q_extra ;
5988 CeedBasisCreateTensorH1Lagrange (ceed , topo_dim , num_comp_u , P , Q ,
60- bp_data .q_mode ,
61- & basis_u );
89+ bp_data .q_mode , & basis_u );
6290 CeedBasisCreateTensorH1Lagrange (ceed , topo_dim , num_comp_x , 2 , Q ,
63- bp_data .q_mode ,
64- & basis_x );
91+ bp_data .q_mode , & basis_x );
6592 CeedBasisGetNumQuadraturePoints (basis_u , & num_qpts );
6693
6794 // CEED restrictions
@@ -155,8 +182,7 @@ PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree,
155182 CeedOperatorSetField (op_setup_rhs , "x" , elem_restr_x , basis_x ,
156183 CEED_VECTOR_ACTIVE );
157184 CeedOperatorSetField (op_setup_rhs , "q_data" , elem_restr_qd_i ,
158- CEED_BASIS_COLLOCATED ,
159- q_data );
185+ CEED_BASIS_COLLOCATED , q_data );
160186 CeedOperatorSetField (op_setup_rhs , "true_soln" , elem_restr_u_i ,
161187 CEED_BASIS_COLLOCATED , * target );
162188 CeedOperatorSetField (op_setup_rhs , "rhs" , elem_restr_u , basis_u ,
@@ -255,3 +281,151 @@ PetscErrorCode CeedLevelTransferSetup(Ceed ceed, CeedInt num_levels,
255281};
256282
257283// -----------------------------------------------------------------------------
284+ // Set up libCEED for BDDC interface vertices
285+ // -----------------------------------------------------------------------------
286+ PetscErrorCode SetupLibceedBDDC (DM dm_vertex , CeedData data_fine ,
287+ CeedDataBDDC data_vertex ,
288+ PetscInt g_vertex_size , PetscInt xl_vertex_size ,
289+ BPData bp_data {
290+ int ierr ;
291+ Ceed ceed = data_fine -> ceed ;
292+ CeedBasis basis_Pi , basis_u = data_fine -> basis_u ;
293+ CeedElemRestriction elem_restr_Pi , elem_restr_Pi_r , elem_restr_r ;
294+ CeedOperator op_Pi_r , op_r_Pi , op_Pi_Pi , op_r_r , op_r_r_inv , op_inject_Pi , op_inject_r , op_restrict_Pi , op_restrict_r ;
295+ CeedVector x_Pi_ceed , y_Pi_ceed , x_r_ceed , y_r_ceed , mask_r_ceed , mask_Gamma_ceed , mask_I_ceed ;
296+ CeedInt topo_dim , num_comp_u , P , Q , num_qpts , num_elem , elem_size ,
297+ q_data_size = bp_data .q_data_size ;
298+
299+ // CEED basis
300+ CeedBasisGetDimension (basis_u , & topo_dim );
301+ CeedBasisGetNumComponents (basis_u , & num_comp_u );
302+ CeedBasisGetNumNodes1D (basis_u , & P );
303+ elem_size = CeedIntPow (P , topo_dim );
304+ CeedBasisGetNumQuadraturePoints1D (basis_u , & Q );
305+ CeedBasisGetNumQuadraturePoints (basis_u , & num_qpts );
306+ CeedScalar * interp_1d , * grad_1d , * q_ref_1d , * q_weight_1d ;
307+ interp_1d = calloc (2 * Q * sizeof (CeedScalar ));
308+ CeedScalar * temp ;
309+ CeedBasisGetInterp1D (basis_u , & temp );
310+ memcpy (interp_1d , temp , Q * sizeof (CeedScalar ));
311+ memcpy (& interp_1d [1 * Q ], temp [(P - 1 )* Q ], Q * sizeof (CeedScalar ));
312+ grad_1d = calloc (2 * Q * sizeof (CeedScalar ));
313+ CeedBasisGetGrad1D (basis_u , & temp );
314+ memcpy (grad_1d , temp , Q * sizeof (CeedScalar ));
315+ memcpy (& grad_1d [1 * Q ], temp [(P - 1 )* Q ], Q * sizeof (CeedScalar ));
316+ q_ref_1d = calloc (Q * sizeof (CeedScalar ));
317+ CeedBasisGetQRef (basis_u , & temp );
318+ memcpy (q_ref_1d , temp , Q * sizeof (CeedScalar ));
319+ q_weight_1d = calloc (Q * sizeof (CeedScalar ));
320+ CeedBasisGetQWeights (basis_u , & temp );
321+ memcpy (q_weight_1d , temp , Q * sizeof (CeedScalar ));
322+ CeedBasisCreateTensorH1 (ceed , topo_dim , num_comp_u , 1 , Q ,
323+ interp_1d , grad_1d , q_ref_1d ,
324+ q_weight_1d , & basis_Pi );
325+
326+ // CEED restrictions
327+ // -- Interface vertex restriction
328+ ierr = CreateRestrictionFromPlex (ceed , dm_vertex , P , topo_dim , 0 , 0 , 0 , & elem_restr_Pi );
329+ CHKERRQ (ierr );
330+
331+ // -- Subdomain restriction
332+ ierr = DMPlexGetHeightStratum (dm_vertex , 0 , & c_start , & c_end ); CHKERRQ (ierr );
333+ num_elem = c_end - c_start ;
334+ CeedInt strides = [num_comp_u , 1 , num_comp_u * elem_size ];
335+ CeedElemRestrictionCreateStrided (ceed , num_elem , elem_size , num_comp_u ,
336+ num_comp_u * num_elem * elem_size , // **NOPAD**
337+ strides , & elem_restr_r );
338+
339+ // -- Subdomain to interface vertex restriction
340+ CeedInt offsets_Pi_r [num_elem * 8 ];
341+ for (CeedInt e = 0 ; e < num_elem ; e ++ ) {
342+ CeedInt elem_stride = e * num_comp_u * elem_size ;
343+ offsets_Pi_r [e * 8 + 0 ] = elem_stride + num_comp_u * 0 ;
344+ offsets_Pi_r [e * 8 + 1 ] = elem_stride + num_comp_u * (P - 1 );
345+ offsets_Pi_r [e * 8 + 2 ] = elem_stride + num_comp_u * (P * (P - 1 ) + 0 );
346+ offsets_Pi_r [e * 8 + 3 ] = elem_stride + num_comp_u * (P * (P - 1 ) + P - 1 );
347+ offsets_Pi_r [e * 8 + 4 ] = elem_stride + num_comp_u * (P * P * (P - 1 ) + 0 );
348+ offsets_Pi_r [e * 8 + 5 ] = elem_stride + num_comp_u * (P * P * (P - 1 ) + P - 1 );
349+ offsets_Pi_r [e * 8 + 6 ] = elem_stride + num_comp_u * (P * P * (P - 1 ) + P * (P - 1 ) + 0 );
350+ offsets_Pi_r [e * 8 + 7 ] = elem_stride + num_comp_u * (P * P * (P - 1 ) + P * (P - 1 ) + P - 1 );
351+ }
352+ CeedElemRestrictionCreate (ceed , num_elem , elem_size , num_comp_u , 1 , num_comp_u * num_elem * elem_size , // **NOPAD**
353+ CEED_MEM_HOST , CEED_COPY_VALUES , & offsets_Pi_r , & elem_restr_Pi_r );
354+
355+ // Create the persistent vectors that will be needed
356+ CeedVectorCreate (ceed , xl_vertex_size , & x_Pi_ceed );
357+ CeedVectorCreate (ceed , xl_vertex_size , & y_Pi_ceed );
358+ CeedVectorCreate (ceed , num_comp_u * elem_size * num_elem , & x_r_ceed );
359+ CeedVectorCreate (ceed , num_comp_u * elem_size * num_elem , & y_r_ceed );
360+
361+ // Create the mass or diff operator
362+ CeedQFunction qf_apply = data_fine -> qf_apply ;
363+ // -- Interface vertices
364+ CeedOperatorCreate (ceed , data_fine -> qf_apply , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_Pi_Pi );
365+ CeedOperatorSetField (op_Pi_Pi , "u" , elem_restr_Pi , basis_u , CEED_VECTOR_ACTIVE );
366+ CeedOperatorSetField (op_Pi_Pi , "q_data" , data_fine -> elem_restr_qd_i ,
367+ CEED_BASIS_COLLOCATED , data_fine -> q_data );
368+ CeedOperatorSetField (op_Pi_Pi , "v" , elem_restr_Pi , basis_u , CEED_VECTOR_ACTIVE );
369+ // -- Subdomains to interface vertices
370+ CeedOperatorCreate (ceed , data_fine -> qf_apply , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_Pi_r );
371+ CeedOperatorSetField (op_Pi_r , "u" , elem_restr_r , basis_u , CEED_VECTOR_ACTIVE );
372+ CeedOperatorSetField (op_Pi_r , "q_data" , data_fine -> elem_restr_qd_i ,
373+ CEED_BASIS_COLLOCATED , data_fine -> q_data );
374+ CeedOperatorSetField (op_Pi_r , "v" , elem_restr_Pi , basis_u , CEED_VECTOR_ACTIVE );
375+ // -- Interface vertices to subdomains
376+ CeedOperatorCreate (ceed , data_fine -> qf_apply , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_r_Pi );
377+ CeedOperatorSetField (op_r_Pi , "u" , elem_restr_Pi , basis_u , CEED_VECTOR_ACTIVE );
378+ CeedOperatorSetField (op_r_Pi , "q_data" , data_fine -> elem_restr_qd_i ,
379+ CEED_BASIS_COLLOCATED , data_fine -> q_data );
380+ CeedOperatorSetField (op_r_Pi , "v" , elem_restr_r , basis_u , CEED_VECTOR_ACTIVE );
381+ // -- Subdomains
382+ CeedOperatorCreate (ceed , data_fine -> qf_apply , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_r_r );
383+ CeedOperatorSetField (op_r_r , "u" , elem_restr_r , basis_u , CEED_VECTOR_ACTIVE );
384+ CeedOperatorSetField (op_r_r , "q_data" , data_fine -> elem_restr_qd_i ,
385+ CEED_BASIS_COLLOCATED , data_fine -> q_data );
386+ CeedOperatorSetField (op_r_r , "v" , elem_restr_r , basis_u , CEED_VECTOR_ACTIVE );
387+ // -- Subdomain FDM inverse
388+ CeedOperatorCreateFDMElementInverse (op_r_r , & op_r_r_inv , CEED_REQUEST_IMMEDIATE );
389+
390+ // Injection and restriction operators
391+ CeedQFunction qf_identity ;
392+ CeedQFunctionCreateIdentity (ceed , num_comp_u , CEED_EVAL_NONE , CEED_EVAL_NONE ,
393+ & qf_identity );
394+ // -- Injection to interface vertices
395+ CeedOperatorCreate (ceed , qf_identity , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_inject_Pi );
396+ CeedOperatorSetField (op_inject_Pi , "input" , elem_restr_r , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE ); // Note: r to Pi
397+ CeedOperatorSetField (op_inject_Pi , "output" , elem_restr_Pi_r , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE );
398+ // -- Injection to subdomains
399+ CeedOperatorCreate (ceed , qf_identity , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_inject_r );
400+ CeedOperatorSetField (op_inject_r , "input" , data_fine -> elem_restr_u , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE );
401+ CeedOperatorSetField (op_inject_r , "output" , elem_restr_r , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE );
402+ // -- Restriction from interface vertices
403+ CeedOperatorCreate (ceed , qf_identity , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_restrict_Pi );
404+ CeedOperatorSetField (op_restrict_Pi , "input" , elem_restr_Pi_r , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE );
405+ CeedOperatorSetField (op_restrict_Pi , "output" , elem_restr_r , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE );
406+ // -- Restriction from subdomains
407+ CeedOperatorCreate (ceed , qf_identity , CEED_QFUNCTION_NONE , CEED_QFUNCTION_NONE , & op_restrict_r );
408+ CeedOperatorSetField (op_restrict_r , "input" , elem_restr_r , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE ); // Note: Pi to r
409+ CeedOperatorSetField (op_restrict_r , "output" , data_fine -> elem_restr_u , CEED_BASIS_COLLOCATED , CEED_VECTOR_ACTIVE );
410+ // -- Cleanup
411+ CeedQFunctionDestroy (& qf_identity );
412+
413+ // Save libCEED data required for level
414+ data_vertex -> basis_Pi = basis_Pi ;
415+ data_vertex -> elem_restr_Pi = elem_restr_Pi ;
416+ data_vertex -> elem_restr_Pi_r = elem_restr_Pi_r ;
417+ data_vertex -> elem_restr_r = elem_restr_r ;
418+ data_vertex -> op_Pi_r = op_Pi_r ;
419+ data_vertex -> op_r_Pi = op_r_Pi ;
420+ data_vertex -> op_Pi_Pi = op_Pi_Pi ;
421+ data_vertex -> op_r_r = op_r_r ;
422+ data_vertex -> op_r_r_inv = op_r_r_inv ;
423+ data_vertex -> x_Pi_ceed = x_Pi_ceed ;
424+ data_vertex -> y_Pi_ceed = y_Pi_ceed ;
425+ data_vertex -> x_r_ceed = x_r_ceed ;
426+ data_vertex -> y_r_ceed = y_r_ceed ;
427+
428+ PetscFunctionReturn (0 );
429+ };
430+
431+ // -----------------------------------------------------------------------------
0 commit comments