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 /// @file 18 /// Diffusion operator example using MFEM 19 #include <ceed.h> 20 #include <mfem.hpp> 21 22 /// A structure used to pass additional data to f_build_diff and f_apply_diff 23 struct BuildContext { CeedInt dim, space_dim; }; 24 25 /// libCEED Q-function for building quadrature data for a diffusion operator 26 static int f_build_diff(void *ctx, CeedInt Q, 27 const CeedScalar *const *in, CeedScalar *const *out) { 28 BuildContext *bc = (BuildContext *)ctx; 29 // in[0] is Jacobians with shape [dim, nc=dim, Q] 30 // in[1] is quadrature weights, size (Q) 31 // 32 // At every quadrature point, compute qw/det(J).adj(J).adj(J)^T and store 33 // the symmetric part of the result. 34 const CeedScalar *J = in[0], *qw = in[1]; 35 CeedScalar *qd = out[0]; 36 switch (bc->dim + 10*bc->space_dim) { 37 case 11: 38 for (CeedInt i=0; i<Q; i++) { 39 qd[i] = qw[i] / J[i]; 40 } 41 break; 42 case 22: 43 for (CeedInt i=0; i<Q; i++) { 44 // J: 0 2 qd: 0 1 adj(J): J22 -J12 45 // 1 3 1 2 -J21 J11 46 const CeedScalar J11 = J[i+Q*0]; 47 const CeedScalar J21 = J[i+Q*1]; 48 const CeedScalar J12 = J[i+Q*2]; 49 const CeedScalar J22 = J[i+Q*3]; 50 const CeedScalar w = qw[i] / (J11*J22 - J21*J12); 51 qd[i+Q*0] = w * (J12*J12 + J22*J22); 52 qd[i+Q*1] = - w * (J11*J12 + J21*J22); 53 qd[i+Q*2] = w * (J11*J11 + J21*J21); 54 } 55 break; 56 case 33: 57 for (CeedInt i=0; i<Q; i++) { 58 // J: 0 3 6 qd: 0 1 2 59 // 1 4 7 1 3 4 60 // 2 5 8 2 4 5 61 const CeedScalar J11 = J[i+Q*0]; 62 const CeedScalar J21 = J[i+Q*1]; 63 const CeedScalar J31 = J[i+Q*2]; 64 const CeedScalar J12 = J[i+Q*3]; 65 const CeedScalar J22 = J[i+Q*4]; 66 const CeedScalar J32 = J[i+Q*5]; 67 const CeedScalar J13 = J[i+Q*6]; 68 const CeedScalar J23 = J[i+Q*7]; 69 const CeedScalar J33 = J[i+Q*8]; 70 const CeedScalar A11 = J22*J33 - J23*J32; 71 const CeedScalar A12 = J13*J32 - J12*J33; 72 const CeedScalar A13 = J12*J23 - J13*J22; 73 const CeedScalar A21 = J23*J31 - J21*J33; 74 const CeedScalar A22 = J11*J33 - J13*J31; 75 const CeedScalar A23 = J13*J21 - J11*J23; 76 const CeedScalar A31 = J21*J32 - J22*J31; 77 const CeedScalar A32 = J12*J31 - J11*J32; 78 const CeedScalar A33 = J11*J22 - J12*J21; 79 const CeedScalar w = qw[i] / (J11*A11 + J21*A12 + J31*A13); 80 qd[i+Q*0] = w * (A11*A11 + A12*A12 + A13*A13); 81 qd[i+Q*1] = w * (A11*A21 + A12*A22 + A13*A23); 82 qd[i+Q*2] = w * (A11*A31 + A12*A32 + A13*A33); 83 qd[i+Q*3] = w * (A21*A21 + A22*A22 + A23*A23); 84 qd[i+Q*4] = w * (A21*A31 + A22*A32 + A23*A33); 85 qd[i+Q*5] = w * (A31*A31 + A32*A32 + A33*A33); 86 } 87 break; 88 default: 89 return CeedError(NULL, 1, "dim=%d, space_dim=%d is not supported", 90 bc->dim, bc->space_dim); 91 } 92 return 0; 93 } 94 95 /// libCEED Q-function for applying a diff operator 96 static int f_apply_diff(void *ctx, CeedInt Q, 97 const CeedScalar *const *in, CeedScalar *const *out) { 98 BuildContext *bc = (BuildContext *)ctx; 99 // in[0], out[0] have shape [dim, nc=1, Q] 100 const CeedScalar *ug = in[0], *qd = in[1]; 101 CeedScalar *vg = out[0]; 102 switch (bc->dim) { 103 case 1: 104 for (CeedInt i=0; i<Q; i++) { 105 vg[i] = ug[i] * qd[i]; 106 } 107 break; 108 case 2: 109 for (CeedInt i=0; i<Q; i++) { 110 const CeedScalar ug0 = ug[i+Q*0]; 111 const CeedScalar ug1 = ug[i+Q*1]; 112 vg[i+Q*0] = qd[i+Q*0]*ug0 + qd[i+Q*1]*ug1; 113 vg[i+Q*1] = qd[i+Q*1]*ug0 + qd[i+Q*2]*ug1; 114 } 115 break; 116 case 3: 117 for (CeedInt i=0; i<Q; i++) { 118 const CeedScalar ug0 = ug[i+Q*0]; 119 const CeedScalar ug1 = ug[i+Q*1]; 120 const CeedScalar ug2 = ug[i+Q*2]; 121 vg[i+Q*0] = qd[i+Q*0]*ug0 + qd[i+Q*1]*ug1 + qd[i+Q*2]*ug2; 122 vg[i+Q*1] = qd[i+Q*1]*ug0 + qd[i+Q*3]*ug1 + qd[i+Q*4]*ug2; 123 vg[i+Q*2] = qd[i+Q*2]*ug0 + qd[i+Q*4]*ug1 + qd[i+Q*5]*ug2; 124 } 125 break; 126 default: 127 return CeedError(NULL, 1, "topo_dim=%d is not supported", bc->dim); 128 } 129 return 0; 130 } 131 132 /// Wrapper for a diffusion CeedOperator as an mfem::Operator 133 class CeedDiffusionOperator : public mfem::Operator { 134 protected: 135 const mfem::FiniteElementSpace *fes; 136 CeedOperator build_oper, oper; 137 CeedBasis basis, mesh_basis; 138 CeedElemRestriction restr, mesh_restr, restr_i, mesh_restr_i; 139 CeedQFunction apply_qfunc, build_qfunc; 140 CeedVector node_coords, rho; 141 142 BuildContext build_ctx; 143 144 CeedVector u, v; 145 146 static void FESpace2Ceed(const mfem::FiniteElementSpace *fes, 147 const mfem::IntegrationRule &ir, 148 Ceed ceed, CeedBasis *basis, 149 CeedElemRestriction *restr) { 150 mfem::Mesh *mesh = fes->GetMesh(); 151 const mfem::FiniteElement *fe = fes->GetFE(0); 152 const int order = fes->GetOrder(0); 153 mfem::Array<int> dof_map; 154 switch (mesh->Dimension()) { 155 case 1: { 156 const mfem::H1_SegmentElement *h1_fe = 157 dynamic_cast<const mfem::H1_SegmentElement *>(fe); 158 MFEM_VERIFY(h1_fe, "invalid FE"); 159 h1_fe->GetDofMap().Copy(dof_map); 160 break; 161 } 162 case 2: { 163 const mfem::H1_QuadrilateralElement *h1_fe = 164 dynamic_cast<const mfem::H1_QuadrilateralElement *>(fe); 165 MFEM_VERIFY(h1_fe, "invalid FE"); 166 h1_fe->GetDofMap().Copy(dof_map); 167 break; 168 } 169 case 3: { 170 const mfem::H1_HexahedronElement *h1_fe = 171 dynamic_cast<const mfem::H1_HexahedronElement *>(fe); 172 MFEM_VERIFY(h1_fe, "invalid FE"); 173 h1_fe->GetDofMap().Copy(dof_map); 174 break; 175 } 176 } 177 const mfem::FiniteElement *fe1d = 178 fes->FEColl()->FiniteElementForGeometry(mfem::Geometry::SEGMENT); 179 mfem::DenseMatrix shape1d(fe1d->GetDof(), ir.GetNPoints()); 180 mfem::DenseMatrix grad1d(fe1d->GetDof(), ir.GetNPoints()); 181 mfem::Vector qref1d(ir.GetNPoints()), qweight1d(ir.GetNPoints()); 182 mfem::Vector shape_i(shape1d.Height()); 183 mfem::DenseMatrix grad_i(grad1d.Height(), 1); 184 const mfem::H1_SegmentElement *h1_fe1d = 185 dynamic_cast<const mfem::H1_SegmentElement *>(fe1d); 186 MFEM_VERIFY(h1_fe1d, "invalid FE"); 187 const mfem::Array<int> &dof_map_1d = h1_fe1d->GetDofMap(); 188 for (int i = 0; i < ir.GetNPoints(); i++) { 189 const mfem::IntegrationPoint &ip = ir.IntPoint(i); 190 qref1d(i) = ip.x; 191 qweight1d(i) = ip.weight; 192 fe1d->CalcShape(ip, shape_i); 193 fe1d->CalcDShape(ip, grad_i); 194 for (int j = 0; j < shape1d.Height(); j++) { 195 shape1d(j,i) = shape_i(dof_map_1d[j]); 196 grad1d(j,i) = grad_i(dof_map_1d[j],0); 197 } 198 } 199 CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes->GetVDim(), order+1, 200 ir.GetNPoints(), shape1d.GetData(), 201 grad1d.GetData(), qref1d.GetData(), 202 qweight1d.GetData(), basis); 203 204 const mfem::Table &el_dof = fes->GetElementToDofTable(); 205 mfem::Array<int> tp_el_dof(el_dof.Size_of_connections()); 206 for (int i = 0; i < mesh->GetNE(); i++) { 207 const int el_offset = fe->GetDof()*i; 208 for (int j = 0; j < fe->GetDof(); j++) { 209 tp_el_dof[j + el_offset] = el_dof.GetJ()[dof_map[j] + el_offset]; 210 } 211 } 212 CeedElemRestrictionCreate(ceed, mesh->GetNE(), fe->GetDof(), 213 fes->GetNDofs(), fes->GetVDim(), CEED_MEM_HOST, CEED_COPY_VALUES, 214 tp_el_dof.GetData(), restr); 215 } 216 217 public: 218 /// Constructor. Assumes @a fes is a scalar FE space. 219 CeedDiffusionOperator(Ceed ceed, const mfem::FiniteElementSpace *fes) 220 : Operator(fes->GetNDofs()), 221 fes(fes) { 222 mfem::Mesh *mesh = fes->GetMesh(); 223 const int order = fes->GetOrder(0); 224 const int ir_order = 2*(order + 2) - 1; // <----- 225 const mfem::IntegrationRule &ir = 226 mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order); 227 CeedInt nqpts, nelem = mesh->GetNE(), dim = mesh->SpaceDimension(); 228 229 FESpace2Ceed(fes, ir, ceed, &basis, &restr); 230 231 const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace(); 232 MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space"); 233 FESpace2Ceed(mesh_fes, ir, ceed, &mesh_basis, &mesh_restr); 234 CeedBasisGetNumQuadraturePoints(basis, &nqpts); 235 236 CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts*dim*(dim+1)/2, 237 nqpts*nelem*dim*(dim+1)/2, 1, &restr_i); 238 CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts, 239 nqpts*nelem, 1, &mesh_restr_i); 240 241 CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &node_coords); 242 CeedVectorSetArray(node_coords, CEED_MEM_HOST, CEED_USE_POINTER, 243 mesh->GetNodes()->GetData()); 244 245 CeedVectorCreate(ceed, nelem*nqpts*dim*(dim+1)/2, &rho); 246 247 // Context data to be passed to the 'f_build_diff' Q-function. 248 build_ctx.dim = mesh->Dimension(); 249 build_ctx.space_dim = mesh->SpaceDimension(); 250 251 // Create the Q-function that builds the diff operator (i.e. computes its 252 // quadrature data) and set its context data. 253 CeedQFunctionCreateInterior(ceed, 1, f_build_diff, 254 __FILE__":f_build_diff", &build_qfunc); 255 CeedQFunctionAddInput(build_qfunc, "dx", dim, CEED_EVAL_GRAD); 256 CeedQFunctionAddInput(build_qfunc, "weights", 1, CEED_EVAL_WEIGHT); 257 CeedQFunctionAddOutput(build_qfunc, "rho", dim*(dim+1)/2, CEED_EVAL_NONE); 258 CeedQFunctionSetContext(build_qfunc, &build_ctx, sizeof(build_ctx)); 259 260 // Create the operator that builds the quadrature data for the diff operator. 261 CeedOperatorCreate(ceed, build_qfunc, NULL, NULL, &build_oper); 262 CeedOperatorSetField(build_oper, "dx", mesh_restr, CEED_NOTRANSPOSE, 263 mesh_basis, CEED_VECTOR_ACTIVE); 264 CeedOperatorSetField(build_oper, "weights", mesh_restr_i, CEED_NOTRANSPOSE, 265 mesh_basis, CEED_VECTOR_NONE); 266 CeedOperatorSetField(build_oper, "rho", restr_i, CEED_NOTRANSPOSE, 267 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 268 269 // Compute the quadrature data for the diff operator. 270 CeedOperatorApply(build_oper, node_coords, rho, 271 CEED_REQUEST_IMMEDIATE); 272 273 // Create the Q-function that defines the action of the diff operator. 274 CeedQFunctionCreateInterior(ceed, 1, f_apply_diff, 275 __FILE__":f_apply_diff", &apply_qfunc); 276 CeedQFunctionAddInput(apply_qfunc, "u", 1, CEED_EVAL_GRAD); 277 CeedQFunctionAddInput(apply_qfunc, "rho", dim*(dim+1)/2, CEED_EVAL_NONE); 278 CeedQFunctionAddOutput(apply_qfunc, "v", 1, CEED_EVAL_GRAD); 279 CeedQFunctionSetContext(apply_qfunc, &build_ctx, sizeof(build_ctx)); 280 281 // Create the diff operator. 282 CeedOperatorCreate(ceed, apply_qfunc, NULL, NULL, &oper); 283 CeedOperatorSetField(oper, "u", restr, CEED_NOTRANSPOSE, 284 basis, CEED_VECTOR_ACTIVE); 285 CeedOperatorSetField(oper, "rho", restr_i, CEED_NOTRANSPOSE, 286 CEED_BASIS_COLLOCATED, rho); 287 CeedOperatorSetField(oper, "v", restr, CEED_NOTRANSPOSE, 288 basis, CEED_VECTOR_ACTIVE); 289 290 CeedVectorCreate(ceed, fes->GetNDofs(), &u); 291 CeedVectorCreate(ceed, fes->GetNDofs(), &v); 292 } 293 294 /// Destructor 295 ~CeedDiffusionOperator() { 296 CeedVectorDestroy(&u); 297 CeedVectorDestroy(&v); 298 CeedVectorDestroy(&rho); 299 CeedVectorDestroy(&node_coords); 300 CeedElemRestrictionDestroy(&restr); 301 CeedElemRestrictionDestroy(&mesh_restr); 302 CeedElemRestrictionDestroy(&restr_i); 303 CeedElemRestrictionDestroy(&mesh_restr_i); 304 CeedBasisDestroy(&basis); 305 CeedBasisDestroy(&mesh_basis); 306 CeedQFunctionDestroy(&build_qfunc); 307 CeedOperatorDestroy(&build_oper); 308 CeedQFunctionDestroy(&apply_qfunc); 309 CeedOperatorDestroy(&oper); 310 } 311 312 /// Operator action 313 virtual void Mult(const mfem::Vector &x, mfem::Vector &y) const { 314 CeedVectorSetArray(u, CEED_MEM_HOST, CEED_USE_POINTER, x.GetData()); 315 CeedVectorSetArray(v, CEED_MEM_HOST, CEED_USE_POINTER, y.GetData()); 316 317 CeedOperatorApply(oper, u, v, CEED_REQUEST_IMMEDIATE); 318 CeedVectorSyncArray(v, CEED_MEM_HOST); 319 } 320 }; 321