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