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 /// libCEED QFunctions for diffusion operator example for a scalar field on the sphere using PETSc 19 20 #ifndef __CUDACC__ 21 # include <math.h> 22 #endif 23 24 // ***************************************************************************** 25 // This QFunction sets up the geometric factors required for integration and 26 // coordinate transformations when reference coordinates have a different 27 // dimension than the one of physical coordinates 28 // 29 // Reference (parent) 2D coordinates: X \in [-1, 1]^2 30 // 31 // Global 3D physical coordinates given by the mesh: xx \in [-R, R]^3 32 // with R radius of the sphere 33 // 34 // Local 3D physical coordinates on the 2D manifold: x \in [-l, l]^3 35 // with l half edge of the cube inscribed in the sphere 36 // 37 // Change of coordinates matrix computed by the library: 38 // (physical 3D coords relative to reference 2D coords) 39 // dxx_j/dX_i (indicial notation) [3 * 2] 40 // 41 // Change of coordinates x (on the 2D manifold) relative to xx (phyisical 3D): 42 // dx_i/dxx_j (indicial notation) [3 * 3] 43 // 44 // Change of coordinates x (on the 2D manifold) relative to X (reference 2D): 45 // (by chain rule) 46 // dx_i/dX_j [3 * 2] = dx_i/dxx_k [3 * 3] * dxx_k/dX_j [3 * 2] 47 // 48 // modJ is given by the magnitude of the cross product of the columns of dx_i/dX_j 49 // 50 // The quadrature data is stored in the array qdata. 51 // 52 // We require the determinant of the Jacobian to properly compute integrals of 53 // the form: int( u v ) 54 // 55 // qdata[0]: modJ * w 56 // 57 // We use the Moore–Penrose (left) pseudoinverse of dx_i/dX_j, to compute dX_i/dx_j (and its transpose), 58 // needed to properly compute integrals of the form: int( gradv gradu ) 59 // 60 // dX_i/dx_j [2 * 3] = (dx_i/dX_j)+ = (dxdX^T dxdX)^(-1) dxdX 61 // 62 // and the product simplifies to yield the contravariant metric tensor 63 // 64 // g^{ij} = dX_i/dx_k dX_j/dx_k = (dxdX^T dxdX)^{-1} 65 // 66 // Stored: g^{ij} (in Voigt convention) in 67 // 68 // qdata[1:3]: [dXdxdXdxT00 dXdxdXdxT01] 69 // [dXdxdXdxT01 dXdxdXdxT11] 70 // ***************************************************************************** 71 72 // ----------------------------------------------------------------------------- 73 CEED_QFUNCTION(SetupDiffGeo)(void *ctx, CeedInt Q, 74 const CeedScalar *const *in, 75 CeedScalar *const *out) { 76 const CeedScalar *X = in[0], *J = in[1], *w = in[2]; 77 CeedScalar *qdata = out[0]; 78 79 // Quadrature Point Loop 80 CeedPragmaSIMD 81 for (CeedInt i=0; i<Q; i++) { 82 // Read global Cartesian coordinates 83 const CeedScalar xx[3] = {X[i+0*Q], 84 X[i+1*Q], 85 X[i+2*Q] 86 }; 87 88 // Read dxxdX Jacobian entries, stored as 89 // 0 3 90 // 1 4 91 // 2 5 92 const CeedScalar dxxdX[3][2] = {{J[i+Q*0], 93 J[i+Q*3]}, 94 {J[i+Q*1], 95 J[i+Q*4]}, 96 {J[i+Q*2], 97 J[i+Q*5]} 98 }; 99 100 // Setup 101 // x = xx (xx^T xx)^{-1/2} 102 // dx/dxx = I (xx^T xx)^{-1/2} - xx xx^T (xx^T xx)^{-3/2} 103 const CeedScalar modxxsq = xx[0]*xx[0]+xx[1]*xx[1]+xx[2]*xx[2]; 104 CeedScalar xxsq[3][3]; 105 for (int j=0; j<3; j++) 106 for (int k=0; k<3; k++) 107 xxsq[j][k] = xx[j]*xx[k] / (sqrt(modxxsq) * modxxsq); 108 109 const CeedScalar dxdxx[3][3] = {{1./sqrt(modxxsq) - xxsq[0][0], 110 -xxsq[0][1], 111 -xxsq[0][2]}, 112 {-xxsq[1][0], 113 1./sqrt(modxxsq) - xxsq[1][1], 114 -xxsq[1][2]}, 115 {-xxsq[2][0], 116 -xxsq[2][1], 117 1./sqrt(modxxsq) - xxsq[2][2]} 118 }; 119 120 CeedScalar dxdX[3][2]; 121 for (int j=0; j<3; j++) 122 for (int k=0; k<2; k++) { 123 dxdX[j][k] = 0; 124 for (int l=0; l<3; l++) 125 dxdX[j][k] += dxdxx[j][l]*dxxdX[l][k]; 126 } 127 128 // J is given by the cross product of the columns of dxdX 129 const CeedScalar J[3]= {dxdX[1][0]*dxdX[2][1] - dxdX[2][0]*dxdX[1][1], 130 dxdX[2][0]*dxdX[0][1] - dxdX[0][0]*dxdX[2][1], 131 dxdX[0][0]*dxdX[1][1] - dxdX[1][0]*dxdX[0][1] 132 }; 133 134 // Use the magnitude of J as our detJ (volume scaling factor) 135 const CeedScalar modJ = sqrt(J[0]*J[0]+J[1]*J[1]+J[2]*J[2]); 136 137 // Interp-to-Interp qdata 138 qdata[i+Q*0] = modJ * w[i]; 139 140 // dxdX_k,j * dxdX_j,k 141 CeedScalar dxdXTdxdX[2][2]; 142 for (int j=0; j<2; j++) 143 for (int k=0; k<2; k++) { 144 dxdXTdxdX[j][k] = 0; 145 for (int l=0; l<3; l++) 146 dxdXTdxdX[j][k] += dxdX[l][j]*dxdX[l][k]; 147 } 148 149 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] 150 -dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 151 152 // Compute inverse of dxdXTdxdX, which is the 2x2 contravariant metric tensor g^{ij} 153 CeedScalar dxdXTdxdXinv[2][2]; 154 dxdXTdxdXinv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 155 dxdXTdxdXinv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 156 dxdXTdxdXinv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 157 dxdXTdxdXinv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 158 159 // Stored in Voigt convention 160 qdata[i+Q*1] = dxdXTdxdXinv[0][0]; 161 qdata[i+Q*2] = dxdXTdxdXinv[1][1]; 162 qdata[i+Q*3] = dxdXTdxdXinv[0][1]; 163 } // End of Quadrature Point Loop 164 165 // Return 166 return 0; 167 } 168 169 // ***************************************************************************** 170 // This QFunction sets up the rhs and true solution for the problem 171 // ***************************************************************************** 172 173 // ----------------------------------------------------------------------------- 174 CEED_QFUNCTION(SetupDiffRhs)(void *ctx, CeedInt Q, 175 const CeedScalar *const *in, 176 CeedScalar *const *out) { 177 // Inputs 178 const CeedScalar *X = in[0], *qdata = in[1]; 179 // Outputs 180 CeedScalar *true_soln = out[0], *rhs = out[1]; 181 182 // Context 183 const CeedScalar *context = (const CeedScalar*)ctx; 184 const CeedScalar R = context[0]; 185 186 // Quadrature Point Loop 187 CeedPragmaSIMD 188 for (CeedInt i=0; i<Q; i++) { 189 // Read global Cartesian coordinates 190 CeedScalar x = X[i+Q*0], y = X[i+Q*1], z = X[i+Q*2]; 191 // Normalize quadrature point coordinates to sphere 192 CeedScalar rad = sqrt(x*x + y*y + z*z); 193 x *= R / rad; 194 y *= R / rad; 195 z *= R / rad; 196 // Compute latitude and longitude 197 const CeedScalar theta = asin(z / R); // latitude 198 const CeedScalar lambda = atan2(y, x); // longitude 199 200 true_soln[i+Q*0] = sin(lambda) * cos(theta); 201 202 rhs[i+Q*0] = qdata[i+Q*0] * 2 * sin(lambda)*cos(theta) / (R*R); 203 204 } // End of Quadrature Point Loop 205 206 return 0; 207 } 208 209 // ***************************************************************************** 210 // This QFunction applies the diffusion operator for a scalar field. 211 // 212 // Inputs: 213 // ug - Input vector gradient at quadrature points 214 // qdata - Geometric factors 215 // 216 // Output: 217 // vg - Output vector (test functions) gradient at quadrature points 218 // 219 // ***************************************************************************** 220 221 // ----------------------------------------------------------------------------- 222 CEED_QFUNCTION(Diff)(void *ctx, CeedInt Q, 223 const CeedScalar *const *in, CeedScalar *const *out) { 224 // Inputs 225 const CeedScalar *ug = in[0], *qdata = in[1]; 226 // Outputs 227 CeedScalar *vg = out[0]; 228 229 // Quadrature Point Loop 230 CeedPragmaSIMD 231 for (CeedInt i=0; i<Q; i++) { 232 // Read spatial derivatives of u 233 const CeedScalar du[2] = {ug[i+Q*0], 234 ug[i+Q*1] 235 }; 236 // Read qdata 237 const CeedScalar wJ = qdata[i+Q*0]; 238 // -- Grad-to-Grad qdata 239 // ---- dXdx_j,k * dXdx_k,j 240 const CeedScalar dXdxdXdxT[2][2] = {{qdata[i+Q*1], 241 qdata[i+Q*3]}, 242 {qdata[i+Q*3], 243 qdata[i+Q*2]} 244 }; 245 246 for (int j=0; j<2; j++) // j = direction of vg 247 vg[i+j*Q] = wJ * (du[0] * dXdxdXdxT[0][j] + 248 du[1] * dXdxdXdxT[1][j]); 249 250 } // End of Quadrature Point Loop 251 252 return 0; 253 } 254 // ----------------------------------------------------------------------------- 255