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 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 #include <ceed.h> 24 #include <math.h> 25 26 // ----------------------------------------------------------------------------- 27 // This QFunction sets up the geometric factors required for integration and 28 // coordinate transformations when reference coordinates have a different 29 // dimension than the one of physical coordinates 30 // 31 // Reference (parent) 2D coordinates: X \in [-1, 1]^2 32 // 33 // Global 3D physical coordinates given by the mesh: xx \in [-R, R]^3 34 // with R radius of the sphere 35 // 36 // Local 3D physical coordinates on the 2D manifold: x \in [-l, l]^3 37 // with l half edge of the cube inscribed in the sphere 38 // 39 // Change of coordinates matrix computed by the library: 40 // (physical 3D coords relative to reference 2D coords) 41 // dxx_j/dX_i (indicial notation) [3 * 2] 42 // 43 // Change of coordinates x (on the 2D manifold) relative to xx (phyisical 3D): 44 // dx_i/dxx_j (indicial notation) [3 * 3] 45 // 46 // Change of coordinates x (on the 2D manifold) relative to X (reference 2D): 47 // (by chain rule) 48 // dx_i/dX_j [3 * 2] = dx_i/dxx_k [3 * 3] * dxx_k/dX_j [3 * 2] 49 // 50 // mod_J is given by the magnitude of the cross product of the columns of dx_i/dX_j 51 // 52 // The quadrature data is stored in the array q_data. 53 // 54 // We require the determinant of the Jacobian to properly compute integrals of 55 // the form: int( u v ) 56 // 57 // q_data[0]: mod_J * w 58 // 59 // We use the Moore–Penrose (left) pseudoinverse of dx_i/dX_j, to compute dX_i/dx_j (and its transpose), 60 // needed to properly compute integrals of the form: int( gradv gradu ) 61 // 62 // dX_i/dx_j [2 * 3] = (dx_i/dX_j)+ = (dxdX^T dxdX)^(-1) dxdX 63 // 64 // and the product simplifies to yield the contravariant metric tensor 65 // 66 // g^{ij} = dX_i/dx_k dX_j/dx_k = (dxdX^T dxdX)^{-1} 67 // 68 // Stored: g^{ij} (in Voigt convention) in 69 // 70 // q_data[1:3]: [dXdxdXdxT00 dXdxdXdxT01] 71 // [dXdxdXdxT01 dXdxdXdxT11] 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 *q_data = 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 mod_xx_sq = xx[0]*xx[0]+xx[1]*xx[1]+xx[2]*xx[2]; 104 CeedScalar xx_sq[3][3]; 105 for (int j=0; j<3; j++) 106 for (int k=0; k<3; k++) 107 xx_sq[j][k] = xx[j]*xx[k] / (sqrt(mod_xx_sq) * mod_xx_sq); 108 109 const CeedScalar dxdxx[3][3] = {{1./sqrt(mod_xx_sq) - xx_sq[0][0], 110 -xx_sq[0][1], 111 -xx_sq[0][2]}, 112 {-xx_sq[1][0], 113 1./sqrt(mod_xx_sq) - xx_sq[1][1], 114 -xx_sq[1][2]}, 115 {-xx_sq[2][0], 116 -xx_sq[2][1], 117 1./sqrt(mod_xx_sq) - xx_sq[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 mod_J = sqrt(J[0]*J[0]+J[1]*J[1]+J[2]*J[2]); 136 137 // Interp-to-Interp q_data 138 q_data[i+Q*0] = mod_J * 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 dxdXTdxdX_inv[2][2]; 154 dxdXTdxdX_inv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 155 dxdXTdxdX_inv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 156 dxdXTdxdX_inv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 157 dxdXTdxdX_inv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 158 159 // Stored in Voigt convention 160 q_data[i+Q*1] = dxdXTdxdX_inv[0][0]; 161 q_data[i+Q*2] = dxdXTdxdX_inv[1][1]; 162 q_data[i+Q*3] = dxdXTdxdX_inv[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 CEED_QFUNCTION(SetupDiffRhs)(void *ctx, CeedInt Q, 173 const CeedScalar *const *in, 174 CeedScalar *const *out) { 175 // Inputs 176 const CeedScalar *X = in[0], *q_data = in[1]; 177 // Outputs 178 CeedScalar *true_soln = out[0], *rhs = out[1]; 179 180 // Context 181 const CeedScalar *context = (const CeedScalar*)ctx; 182 const CeedScalar R = context[0]; 183 184 // Quadrature Point Loop 185 CeedPragmaSIMD 186 for (CeedInt i=0; i<Q; i++) { 187 // Read global Cartesian coordinates 188 CeedScalar x = X[i+Q*0], y = X[i+Q*1], z = X[i+Q*2]; 189 // Normalize quadrature point coordinates to sphere 190 CeedScalar rad = sqrt(x*x + y*y + z*z); 191 x *= R / rad; 192 y *= R / rad; 193 z *= R / rad; 194 // Compute latitude and longitude 195 const CeedScalar theta = asin(z / R); // latitude 196 const CeedScalar lambda = atan2(y, x); // longitude 197 198 true_soln[i+Q*0] = sin(lambda) * cos(theta); 199 200 rhs[i+Q*0] = q_data[i+Q*0] * 2 * sin(lambda)*cos(theta) / (R*R); 201 202 } // End of Quadrature Point Loop 203 204 return 0; 205 } 206 207 // ----------------------------------------------------------------------------- 208 // This QFunction applies the diffusion operator for a scalar field. 209 // 210 // Inputs: 211 // ug - Input vector gradient at quadrature points 212 // q_data - Geometric factors 213 // 214 // Output: 215 // vg - Output vector (test functions) gradient at quadrature points 216 // 217 // ----------------------------------------------------------------------------- 218 CEED_QFUNCTION(Diff)(void *ctx, CeedInt Q, 219 const CeedScalar *const *in, CeedScalar *const *out) { 220 // Inputs 221 const CeedScalar *ug = in[0], *q_data = in[1]; 222 // Outputs 223 CeedScalar *vg = out[0]; 224 225 // Quadrature Point Loop 226 CeedPragmaSIMD 227 for (CeedInt i=0; i<Q; i++) { 228 // Read spatial derivatives of u 229 const CeedScalar du[2] = {ug[i+Q*0], 230 ug[i+Q*1] 231 }; 232 // Read q_data 233 const CeedScalar w_det_J = q_data[i+Q*0]; 234 // -- Grad-to-Grad q_data 235 // ---- dXdx_j,k * dXdx_k,j 236 const CeedScalar dXdxdXdx_T[2][2] = {{q_data[i+Q*1], 237 q_data[i+Q*3]}, 238 {q_data[i+Q*3], 239 q_data[i+Q*2]} 240 }; 241 242 for (int j=0; j<2; j++) // j = direction of vg 243 vg[i+j*Q] = w_det_J * (du[0] * dXdxdXdx_T[0][j] + 244 du[1] * dXdxdXdx_T[1][j]); 245 246 } // End of Quadrature Point Loop 247 248 return 0; 249 } 250 // ----------------------------------------------------------------------------- 251 252 #endif // bp3sphere_h 253