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 // Stored: dX_i/dx_j * dX_j/dx_i (in Voigt convention) 63 // in qdata[1:3] as 64 // [dXdxdXdxT11 dXdxdXdxT12] 65 // [dXdxdXdxT21 dXdxdXdxT22] 66 // ***************************************************************************** 67 68 // ----------------------------------------------------------------------------- 69 CEED_QFUNCTION(SetupDiffGeo)(void *ctx, CeedInt Q, 70 const CeedScalar *const *in, 71 CeedScalar *const *out) { 72 const CeedScalar *X = in[0], *J = in[1], *w = in[2]; 73 CeedScalar *qdata = out[0]; 74 75 // Quadrature Point Loop 76 CeedPragmaSIMD 77 for (CeedInt i=0; i<Q; i++) { 78 // Read global Cartesian coordinates 79 const CeedScalar xx[3] = {X[i+0*Q], 80 X[i+1*Q], 81 X[i+2*Q] 82 }; 83 84 // Read dxxdX Jacobian entries, stored as 85 // 0 3 86 // 1 4 87 // 2 5 88 const CeedScalar dxxdX[3][2] = {{J[i+Q*0], 89 J[i+Q*3]}, 90 {J[i+Q*1], 91 J[i+Q*4]}, 92 {J[i+Q*2], 93 J[i+Q*5]} 94 }; 95 96 // Setup 97 // x = xx (xx^T xx)^{-1/2} 98 // dx/dxx = I (xx^T xx)^{-1/2} - xx xx^T (xx^T xx)^{-3/2} 99 const CeedScalar modxxsq = xx[0]*xx[0]+xx[1]*xx[1]+xx[2]*xx[2]; 100 CeedScalar xxsq[3][3]; 101 for (int j=0; j<3; j++) 102 for (int k=0; k<3; k++) 103 xxsq[j][k] = xx[j]*xx[k] / (sqrt(modxxsq) * modxxsq); 104 105 const CeedScalar dxdxx[3][3] = {{1./sqrt(modxxsq) - xxsq[0][0], 106 -xxsq[0][1], 107 -xxsq[0][2]}, 108 {-xxsq[1][0], 109 1./sqrt(modxxsq) - xxsq[1][1], 110 -xxsq[1][2]}, 111 {-xxsq[2][0], 112 -xxsq[2][1], 113 1./sqrt(modxxsq) - xxsq[2][2]} 114 }; 115 116 CeedScalar dxdX[3][2]; 117 for (int j=0; j<3; j++) 118 for (int k=0; k<2; k++) { 119 dxdX[j][k] = 0; 120 for (int l=0; l<3; l++) 121 dxdX[j][k] += dxdxx[j][l]*dxxdX[l][k]; 122 } 123 124 // J is given by the cross product of the columns of dxdX 125 const CeedScalar J[3]= {dxdX[1][0]*dxdX[2][1] - dxdX[2][0]*dxdX[1][1], 126 dxdX[2][0]*dxdX[0][1] - dxdX[0][0]*dxdX[2][1], 127 dxdX[0][0]*dxdX[1][1] - dxdX[1][0]*dxdX[0][1] 128 }; 129 130 // Use the magnitude of J as our detJ (volume scaling factor) 131 const CeedScalar modJ = sqrt(J[0]*J[0]+J[1]*J[1]+J[2]*J[2]); 132 133 // Interp-to-Interp qdata 134 qdata[i+Q*0] = modJ * w[i]; 135 136 // dxdX_j,k * dxdX_k,j, needed for the pseudoinverse 137 CeedScalar dxdXTdxdX[2][2]; 138 for (int j=0; j<2; j++) 139 for (int k=0; k<2; k++) { 140 dxdXTdxdX[j][k] = 0; 141 for (int l=0; l<3; l++) 142 dxdXTdxdX[j][k] += dxdX[l][j]*dxdX[l][k]; 143 } 144 145 const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1] 146 -dxdXTdxdX[1][0] * dxdXTdxdX[0][1]; 147 148 // Compute inverse of dxdXTdxdX, needed for the pseudoinverse 149 CeedScalar dxdXTdxdXinv[2][2]; 150 dxdXTdxdXinv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX; 151 dxdXTdxdXinv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX; 152 dxdXTdxdXinv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX; 153 dxdXTdxdXinv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX; 154 155 // Compute the pseudo inverse of dxdX 156 CeedScalar pseudodXdx[2][3]; 157 for (int j=0; j<2; j++) 158 for (int k=0; k<3; k++) { 159 pseudodXdx[j][k] = 0; 160 for (int l=0; l<2; l++) 161 pseudodXdx[j][k] += dxdXTdxdXinv[j][l]*dxdX[k][l]; 162 } 163 164 // Grad-to-Grad qdata is given by pseudodXdx * pseudodXdxT 165 CeedScalar dXdxdXdxT[2][2]; 166 for (int j=0; j<2; j++) 167 for (int k=0; k<2; k++) { 168 dXdxdXdxT[j][k] = 0; 169 for (int l=0; l<3; l++) 170 dXdxdXdxT[j][k] += pseudodXdx[j][l]*pseudodXdx[k][l]; 171 } 172 173 // Stored in Voigt convention 174 qdata[i+Q*1] = dXdxdXdxT[0][0]; 175 qdata[i+Q*2] = dXdxdXdxT[1][1]; 176 qdata[i+Q*3] = dXdxdXdxT[0][1]; 177 178 } // End of Quadrature Point Loop 179 180 // Return 181 return 0; 182 } 183 184 // ***************************************************************************** 185 // This QFunction sets up the rhs and true solution for the problem 186 // ***************************************************************************** 187 188 // ----------------------------------------------------------------------------- 189 CEED_QFUNCTION(SetupDiffRhs)(void *ctx, CeedInt Q, 190 const CeedScalar *const *in, 191 CeedScalar *const *out) { 192 // Inputs 193 const CeedScalar *X = in[0], *qdata = in[1]; 194 // Outputs 195 CeedScalar *true_soln = out[0], *rhs = out[1]; 196 197 // Context 198 const CeedScalar *context = (const CeedScalar*)ctx; 199 const CeedScalar R = context[0]; 200 201 // Quadrature Point Loop 202 CeedPragmaSIMD 203 for (CeedInt i=0; i<Q; i++) { 204 // Read global Cartesian coordinates 205 CeedScalar x = X[i+Q*0], y = X[i+Q*1], z = X[i+Q*2]; 206 // Normalize quadrature point coordinates to sphere 207 CeedScalar rad = sqrt(x*x + y*y + z*z); 208 x *= R / rad; 209 y *= R / rad; 210 z *= R / rad; 211 // Compute latitude and longitude 212 const CeedScalar theta = asin(z / R); // latitude 213 const CeedScalar lambda = atan2(y, x); // longitude 214 215 true_soln[i+Q*0] = sin(lambda) * cos(theta); 216 217 rhs[i+Q*0] = qdata[i+Q*0] * 2 * sin(lambda)*cos(theta) / (R*R); 218 219 } // End of Quadrature Point Loop 220 221 return 0; 222 } 223 224 // ***************************************************************************** 225 // This QFunction applies the diffusion operator for a scalar field. 226 // 227 // Inputs: 228 // ug - Input vector gradient at quadrature points 229 // qdata - Geometric factors 230 // 231 // Output: 232 // vg - Output vector (test functions) gradient at quadrature points 233 // 234 // ***************************************************************************** 235 236 // ----------------------------------------------------------------------------- 237 CEED_QFUNCTION(Diff)(void *ctx, CeedInt Q, 238 const CeedScalar *const *in, CeedScalar *const *out) { 239 // Inputs 240 const CeedScalar *ug = in[0], *qdata = in[1]; 241 // Outputs 242 CeedScalar *vg = out[0]; 243 244 // Quadrature Point Loop 245 CeedPragmaSIMD 246 for (CeedInt i=0; i<Q; i++) { 247 // Read spatial derivatives of u 248 const CeedScalar du[2] = {ug[i+Q*0], 249 ug[i+Q*1] 250 }; 251 // Read qdata 252 const CeedScalar wJ = qdata[i+Q*0]; 253 // -- Grad-to-Grad qdata 254 // ---- dXdx_j,k * dXdx_k,j 255 const CeedScalar dXdxdXdxT[2][2] = {{qdata[i+Q*1], 256 qdata[i+Q*3]}, 257 {qdata[i+Q*3], 258 qdata[i+Q*2]} 259 }; 260 261 for (int j=0; j<2; j++) // j = direction of vg 262 vg[i+j*Q] = wJ * (du[0] * dXdxdXdxT[0][j] + 263 du[1] * dXdxdXdxT[1][j]); 264 265 } // End of Quadrature Point Loop 266 267 return 0; 268 } 269 // ----------------------------------------------------------------------------- 270