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