1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include "../include/libceedsetup.h" 9 10 #include <stdio.h> 11 12 #include "../include/petscutils.h" 13 14 // ----------------------------------------------------------------------------- 15 // Destroy libCEED operator objects 16 // ----------------------------------------------------------------------------- 17 PetscErrorCode CeedDataDestroy(CeedInt i, CeedData data) { 18 PetscFunctionBeginUser; 19 CeedVectorDestroy(&data->q_data); 20 CeedVectorDestroy(&data->x_ceed); 21 CeedVectorDestroy(&data->y_ceed); 22 CeedBasisDestroy(&data->basis_x); 23 CeedBasisDestroy(&data->basis_u); 24 CeedElemRestrictionDestroy(&data->elem_restr_u); 25 CeedElemRestrictionDestroy(&data->elem_restr_x); 26 CeedElemRestrictionDestroy(&data->elem_restr_u_i); 27 CeedElemRestrictionDestroy(&data->elem_restr_qd_i); 28 CeedQFunctionDestroy(&data->qf_apply); 29 CeedOperatorDestroy(&data->op_apply); 30 if (i > 0) { 31 CeedOperatorDestroy(&data->op_prolong); 32 CeedOperatorDestroy(&data->op_restrict); 33 } 34 PetscCall(PetscFree(data)); 35 PetscFunctionReturn(PETSC_SUCCESS); 36 }; 37 38 // ----------------------------------------------------------------------------- 39 // Set up libCEED for a given degree 40 // ----------------------------------------------------------------------------- 41 PetscErrorCode SetupLibceedByDegree(DM dm, Ceed ceed, CeedInt degree, CeedInt topo_dim, CeedInt q_extra, PetscInt num_comp_x, PetscInt num_comp_u, 42 PetscInt g_size, PetscInt xl_size, BPData bp_data, CeedData data, PetscBool setup_rhs, CeedVector rhs_ceed, 43 CeedVector *target) { 44 DM dm_coord; 45 Vec coords; 46 const PetscScalar *coord_array; 47 CeedBasis basis_x, basis_u; 48 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i; 49 CeedQFunction qf_setup_geo, qf_apply; 50 CeedOperator op_setup_geo, op_apply; 51 CeedVector x_coord, q_data, x_ceed, y_ceed; 52 CeedInt num_qpts, c_start, c_end, num_elem, q_data_size = bp_data.q_data_size; 53 CeedScalar R = 1; // radius of the sphere 54 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube 55 56 PetscFunctionBeginUser; 57 PetscCall(DMGetCoordinateDM(dm, &dm_coord)); 58 59 // CEED bases 60 PetscCall(CreateBasisFromPlex(ceed, dm_coord, 0, 0, 0, 0, bp_data, &basis_x)); 61 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 62 63 // CEED restrictions 64 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, 0, 0, &elem_restr_x)); 65 PetscCall(CreateRestrictionFromPlex(ceed, dm, 0, 0, 0, &elem_restr_u)); 66 67 PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); 68 num_elem = c_end - c_start; 69 CeedBasisGetNumQuadraturePoints(basis_u, &num_qpts); 70 71 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, num_comp_u * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_u_i); 72 CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size, q_data_size * num_elem * num_qpts, CEED_STRIDES_BACKEND, &elem_restr_qd_i); 73 74 // Element coordinates 75 PetscCall(DMGetCoordinatesLocal(dm, &coords)); 76 PetscCall(VecGetArrayRead(coords, &coord_array)); 77 78 CeedElemRestrictionCreateVector(elem_restr_x, &x_coord, NULL); 79 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (PetscScalar *)coord_array); 80 PetscCall(VecRestoreArrayRead(coords, &coord_array)); 81 82 // Create the persistent vectors that will be needed in setup and apply 83 CeedVectorCreate(ceed, q_data_size * num_elem * num_qpts, &q_data); 84 CeedVectorCreate(ceed, xl_size, &x_ceed); 85 CeedVectorCreate(ceed, xl_size, &y_ceed); 86 87 // Create the QFunction that builds the context data 88 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup_geo); 89 CeedQFunctionAddInput(qf_setup_geo, "x", num_comp_x, CEED_EVAL_INTERP); 90 CeedQFunctionAddInput(qf_setup_geo, "dx", num_comp_x * topo_dim, CEED_EVAL_GRAD); 91 CeedQFunctionAddInput(qf_setup_geo, "weight", 1, CEED_EVAL_WEIGHT); 92 CeedQFunctionAddOutput(qf_setup_geo, "qdata", q_data_size, CEED_EVAL_NONE); 93 94 // Create the operator that builds the quadrature data 95 CeedOperatorCreate(ceed, qf_setup_geo, NULL, NULL, &op_setup_geo); 96 CeedOperatorSetField(op_setup_geo, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 97 CeedOperatorSetField(op_setup_geo, "dx", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 98 CeedOperatorSetField(op_setup_geo, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 99 CeedOperatorSetField(op_setup_geo, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE); 100 101 // Setup q_data 102 CeedOperatorApply(op_setup_geo, x_coord, q_data, CEED_REQUEST_IMMEDIATE); 103 104 // Set up PDE operator 105 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? topo_dim : 1; 106 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? topo_dim : 1; 107 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply); 108 CeedQFunctionAddInput(qf_apply, "u", num_comp_u * in_scale, bp_data.in_mode); 109 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE); 110 CeedQFunctionAddOutput(qf_apply, "v", num_comp_u * out_scale, bp_data.out_mode); 111 112 // Create the mass or diff operator 113 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 114 CeedOperatorSetField(op_apply, "u", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 115 CeedOperatorSetField(op_apply, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 116 CeedOperatorSetField(op_apply, "v", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 117 118 // Set up RHS if needed 119 if (setup_rhs) { 120 CeedQFunction qf_setup_rhs; 121 CeedOperator op_setup_rhs; 122 CeedVectorCreate(ceed, num_elem * num_qpts * num_comp_u, target); 123 // Create the q-function that sets up the RHS and true solution 124 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs); 125 CeedQFunctionAddInput(qf_setup_rhs, "x", num_comp_x, CEED_EVAL_INTERP); 126 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE); 127 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp_u, CEED_EVAL_NONE); 128 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp_u, CEED_EVAL_INTERP); 129 130 // Create the operator that builds the RHS and true solution 131 CeedOperatorCreate(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs); 132 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); 133 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_qd_i, CEED_BASIS_NONE, q_data); 134 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_i, CEED_BASIS_NONE, *target); 135 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u, basis_u, CEED_VECTOR_ACTIVE); 136 137 // Set up the libCEED context 138 CeedQFunctionContext ctx_rhs_setup; 139 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup); 140 CeedScalar rhs_setup_data[2] = {R, l}; 141 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data); 142 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup); 143 CeedQFunctionContextDestroy(&ctx_rhs_setup); 144 145 // Setup RHS and target 146 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE); 147 148 // Cleanup 149 CeedQFunctionDestroy(&qf_setup_rhs); 150 CeedOperatorDestroy(&op_setup_rhs); 151 } 152 153 // Cleanup 154 CeedQFunctionDestroy(&qf_setup_geo); 155 CeedOperatorDestroy(&op_setup_geo); 156 CeedVectorDestroy(&x_coord); 157 158 // Save libCEED data required for level 159 data->basis_x = basis_x; 160 data->basis_u = basis_u; 161 data->elem_restr_x = elem_restr_x; 162 data->elem_restr_u = elem_restr_u; 163 data->elem_restr_u_i = elem_restr_u_i; 164 data->elem_restr_qd_i = elem_restr_qd_i; 165 data->qf_apply = qf_apply; 166 data->op_apply = op_apply; 167 data->q_data = q_data; 168 data->x_ceed = x_ceed; 169 data->y_ceed = y_ceed; 170 data->q_data_size = q_data_size; 171 PetscFunctionReturn(PETSC_SUCCESS); 172 }; 173 174 // ----------------------------------------------------------------------------- 175 // Setup libCEED level transfer operator objects 176 // ----------------------------------------------------------------------------- 177 PetscErrorCode CeedLevelTransferSetup(DM dm, Ceed ceed, CeedInt level, CeedInt num_comp_u, CeedData *data, BPData bp_data, Vec fine_mult) { 178 PetscFunctionBeginUser; 179 // Restriction - Fine to corse 180 CeedOperator op_restrict; 181 // Interpolation - Corse to fine 182 CeedOperator op_prolong; 183 // Coarse grid operator 184 CeedOperator op_apply; 185 // Basis 186 CeedBasis basis_u; 187 PetscCall(CreateBasisFromPlex(ceed, dm, 0, 0, 0, 0, bp_data, &basis_u)); 188 189 // --------------------------------------------------------------------------- 190 // Coarse Grid, Prolongation, and Restriction Operators 191 // --------------------------------------------------------------------------- 192 // Create the Operators that compute the prolongation and 193 // restriction between the p-multigrid levels and the coarse grid eval. 194 // --------------------------------------------------------------------------- 195 // Place in libCEED array 196 PetscMemType m_mem_type; 197 PetscCall(VecReadP2C(fine_mult, &m_mem_type, data[level]->x_ceed)); 198 199 CeedOperatorMultigridLevelCreate(data[level]->op_apply, data[level]->x_ceed, data[level - 1]->elem_restr_u, basis_u, &op_apply, &op_prolong, 200 &op_restrict); 201 202 // Restore PETSc vector 203 PetscCall(VecReadC2P(data[level]->x_ceed, m_mem_type, fine_mult)); 204 PetscCall(VecZeroEntries(fine_mult)); 205 // -- Save libCEED data 206 data[level - 1]->op_apply = op_apply; 207 data[level]->op_prolong = op_prolong; 208 data[level]->op_restrict = op_restrict; 209 210 CeedBasisDestroy(&basis_u); 211 PetscFunctionReturn(PETSC_SUCCESS); 212 }; 213 214 // ----------------------------------------------------------------------------- 215