1 static char help[] = "Biological network from https://link.springer.com/article/10.1007/s42967-023-00297-3\n\n\n"; 2 3 #include <petscts.h> 4 #include <petscsf.h> 5 #include <petscdmplex.h> 6 #include <petscdmplextransform.h> 7 #include <petscdmforest.h> 8 #include <petscviewerhdf5.h> 9 #include <petscds.h> 10 11 /* 12 Here we solve the system of PDEs on \Omega \in R^d: 13 14 * dC/dt - D^2 \Delta C - \nabla p \cross \nabla p + \alpha sqrt(||C||^2_F + eps)^(\gamma-2) C = 0 15 * - \nabla \cdot ((r + C) \nabla p) = S 16 17 where: 18 C = symmetric dxd conductivity tensor 19 p = potential 20 S = source 21 22 with natural boundary conditions on \partial\Omega: 23 \nabla C \cdot n = 0 24 \nabla ((r + C)\nabla p) \cdot n = 0 25 26 Parameters: 27 D = diffusion constant 28 \alpha = metabolic coefficient 29 \gamma = metabolic exponent 30 r, eps are regularization parameters 31 32 We use Lagrange elements for C_ij and P. 33 Equations are rescaled to obtain a symmetric Jacobian. 34 */ 35 36 typedef enum _fieldidx { 37 C_FIELD_ID = 0, 38 P_FIELD_ID, 39 NUM_FIELDS 40 } FieldIdx; 41 42 typedef enum _constantidx { 43 R_ID = 0, 44 EPS_ID, 45 ALPHA_ID, 46 GAMMA_ID, 47 D_ID, 48 FIXC_ID, 49 SPLIT_ID, 50 NUM_CONSTANTS 51 } ConstantIdx; 52 53 PetscLogStage SetupStage, SolveStage; 54 55 #define NORM2C(c00, c01, c11) (PetscSqr(c00) + 2 * PetscSqr(c01) + PetscSqr(c11)) 56 #define NORM2C_3d(c00, c01, c02, c11, c12, c22) (PetscSqr(c00) + 2 * PetscSqr(c01) + 2 * PetscSqr(c02) + PetscSqr(c11) + 2 * PetscSqr(c12) + PetscSqr(c22)) 57 58 /* Eigenvalues real 3x3 symmetric matrix https://onlinelibrary.wiley.com/doi/full/10.1002/nme.7153 */ 59 #if !PetscDefined(USE_COMPLEX) 60 static inline void Eigenvalues_Sym3x3(PetscReal a11, PetscReal a12, PetscReal a13, PetscReal a22, PetscReal a23, PetscReal a33, PetscReal x[3]) 61 { 62 const PetscBool td = (PetscBool)(a13 == 0 && a23 == 0); 63 if (td) { /* two-dimensional case */ 64 PetscReal a12_2 = PetscSqr(a12); 65 PetscReal delta = PetscSqr(a11 - a22) + 4 * a12_2; 66 PetscReal b = -(a11 + a22); 67 PetscReal c = a11 * a22 - a12_2; 68 PetscReal temp = -0.5 * (b + PetscCopysignReal(1.0, b) * PetscSqrtReal(delta)); 69 70 x[0] = temp; 71 x[1] = c / temp; 72 x[2] = a33; 73 } else { 74 const PetscReal I1 = a11 + a22 + a33; 75 const PetscReal J2 = (PetscSqr(a11 - a22) + PetscSqr(a22 - a33) + PetscSqr(a33 - a11)) / 6 + PetscSqr(a12) + PetscSqr(a23) + PetscSqr(a13); 76 const PetscReal s = PetscSqrtReal(J2 / 3); 77 const PetscReal tol = PETSC_MACHINE_EPSILON; 78 79 if (s < tol) { 80 x[0] = x[1] = x[2] = 0.0; 81 } else { 82 const PetscReal S[] = {a11 - I1 / 3, a12, a13, a22 - I1 / 3, a23, a33 - I1 / 3}; 83 84 /* T = S^2 */ 85 PetscReal T[6]; 86 T[0] = S[0] * S[0] + S[1] * S[1] + S[2] * S[2]; 87 T[1] = S[0] * S[1] + S[1] * S[3] + S[2] * S[4]; 88 T[2] = S[0] * S[2] + S[1] * S[4] + S[2] * S[5]; 89 T[3] = S[1] * S[1] + S[3] * S[3] + S[4] * S[4]; 90 T[4] = S[1] * S[2] + S[3] * S[4] + S[4] * S[5]; 91 T[5] = S[2] * S[2] + S[4] * S[4] + S[5] * S[5]; 92 93 T[0] = T[0] - 2.0 * J2 / 3.0; 94 T[3] = T[3] - 2.0 * J2 / 3.0; 95 T[5] = T[5] - 2.0 * J2 / 3.0; 96 97 const PetscReal aa = NORM2C_3d(T[0] - s * S[0], T[1] - s * S[1], T[2] - s * S[2], T[3] - s * S[3], T[4] - s * S[4], T[5] - s * S[5]); 98 const PetscReal bb = NORM2C_3d(T[0] + s * S[0], T[1] + s * S[1], T[2] + s * S[2], T[3] + s * S[3], T[4] + s * S[4], T[5] + s * S[5]); 99 const PetscReal d = PetscSqrtReal(aa / bb); 100 const PetscBool sj = (PetscBool)(1.0 - d > 0.0); 101 102 if (PetscAbsReal(1 - d) < tol) { 103 x[0] = -PetscSqrtReal(J2); 104 x[1] = 0.0; 105 x[2] = PetscSqrtReal(J2); 106 } else { 107 const PetscReal sjn = sj ? 1.0 : -1.0; 108 //const PetscReal atanarg = sj ? d : 1.0 / d; 109 //const PetscReal alpha = 2.0 * PetscAtanReal(atanarg) / 3.0; 110 const PetscReal atanval = d > tol ? 2.0 * PetscAtanReal(sj ? d : 1.0 / d) : (sj ? 0.0 : PETSC_PI); 111 const PetscReal alpha = atanval / 3.0; 112 const PetscReal cd = s * PetscCosReal(alpha) * sjn; 113 const PetscReal sd = PetscSqrtReal(J2) * PetscSinReal(alpha); 114 115 x[0] = 2 * cd; 116 x[1] = -cd + sd; 117 x[2] = -cd - sd; 118 } 119 } 120 x[0] += I1 / 3.0; 121 x[1] += I1 / 3.0; 122 x[2] += I1 / 3.0; 123 } 124 } 125 #endif 126 127 /* compute shift to make C positive definite */ 128 static inline PetscReal FIX_C_3d(PetscScalar C00, PetscScalar C01, PetscScalar C02, PetscScalar C11, PetscScalar C12, PetscScalar C22) 129 { 130 #if !PetscDefined(USE_COMPLEX) 131 PetscReal eigs[3], s = 0.0; 132 PetscBool twod = (PetscBool)(C02 == 0 && C12 == 0 && C22 == 0); 133 Eigenvalues_Sym3x3(C00, C01, C02, C11, C12, C22, eigs); 134 if (twod) eigs[2] = 1.0; 135 if (eigs[0] <= 0 || eigs[1] <= 0 || eigs[2] <= 0) s = -PetscMin(eigs[0], PetscMin(eigs[1], eigs[2])) + PETSC_SMALL; 136 return s; 137 #else 138 return 0.0; 139 #endif 140 } 141 142 static inline PetscReal FIX_C(PetscScalar C00, PetscScalar C01, PetscScalar C11) 143 { 144 return FIX_C_3d(C00, C01, 0, C11, 0, 0); 145 } 146 147 /* residual for C when tested against basis functions */ 148 static void C_0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 149 { 150 const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]); 151 const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]); 152 const PetscReal eps = PetscRealPart(constants[EPS_ID]); 153 const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0); 154 const PetscScalar *gradp = split ? a_x + aOff_x[P_FIELD_ID] : u_x + uOff_x[P_FIELD_ID]; 155 const PetscScalar outerp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[1] * gradp[1]}; 156 const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]]; 157 const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1]; 158 const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2]; 159 const PetscScalar norm = NORM2C(C00, C01, C11) + eps; 160 const PetscScalar nexp = (gamma - 2.0) / 2.0; 161 const PetscScalar fnorm = PetscPowScalar(norm, nexp); 162 const PetscScalar eqss[] = {0.5, 1., 0.5}; /* equations rescaling for a symmetric Jacobian */ 163 164 for (PetscInt k = 0; k < 3; k++) f0[k] = eqss[k] * (u_t[uOff[C_FIELD_ID] + k] - outerp[k] + alpha * fnorm * u[uOff[C_FIELD_ID] + k]); 165 } 166 167 /* 3D version */ 168 static void C_0_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 169 { 170 const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]); 171 const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]); 172 const PetscReal eps = PetscRealPart(constants[EPS_ID]); 173 const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0); 174 const PetscScalar *gradp = split ? a_x + aOff_x[P_FIELD_ID] : u_x + uOff_x[P_FIELD_ID]; 175 const PetscScalar outerp[] = {gradp[0] * gradp[0], gradp[0] * gradp[1], gradp[0] * gradp[2], gradp[1] * gradp[1], gradp[1] * gradp[2], gradp[2] * gradp[2]}; 176 const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]]; 177 const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1]; 178 const PetscScalar C02 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2]; 179 const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 3] : u[uOff[C_FIELD_ID] + 3]; 180 const PetscScalar C12 = split ? a[aOff[C_FIELD_ID] + 4] : u[uOff[C_FIELD_ID] + 4]; 181 const PetscScalar C22 = split ? a[aOff[C_FIELD_ID] + 5] : u[uOff[C_FIELD_ID] + 5]; 182 const PetscScalar norm = NORM2C_3d(C00, C01, C02, C11, C12, C22) + eps; 183 const PetscScalar nexp = (gamma - 2.0) / 2.0; 184 const PetscScalar fnorm = PetscPowScalar(norm, nexp); 185 const PetscScalar eqss[] = {0.5, 1., 1., 0.5, 1., 0.5}; 186 187 for (PetscInt k = 0; k < 6; k++) f0[k] = eqss[k] * (u_t[uOff[C_FIELD_ID] + k] - outerp[k] + alpha * fnorm * u[uOff[C_FIELD_ID] + k]); 188 } 189 190 /* Jacobian for C against C basis functions */ 191 static void JC_0_c0c0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 192 { 193 const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]); 194 const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]); 195 const PetscReal eps = PetscRealPart(constants[EPS_ID]); 196 const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0); 197 const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]]; 198 const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1]; 199 const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2]; 200 const PetscScalar norm = NORM2C(C00, C01, C11) + eps; 201 const PetscScalar nexp = (gamma - 2.0) / 2.0; 202 const PetscScalar fnorm = PetscPowScalar(norm, nexp); 203 const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0); 204 const PetscScalar dC[] = {2 * C00, 4 * C01, 2 * C11}; 205 const PetscScalar eqss[] = {0.5, 1., 0.5}; 206 207 for (PetscInt k = 0; k < 3; k++) { 208 if (!split) { 209 for (PetscInt j = 0; j < 3; j++) J[k * 3 + j] = eqss[k] * (alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k]); 210 } 211 J[k * 3 + k] += eqss[k] * (alpha * fnorm + u_tShift); 212 } 213 } 214 215 static void JC_0_c0c0_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 216 { 217 const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]); 218 const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]); 219 const PetscReal eps = PetscRealPart(constants[EPS_ID]); 220 const PetscBool split = (PetscBool)(PetscRealPart(constants[SPLIT_ID]) != 0.0); 221 const PetscScalar C00 = split ? a[aOff[C_FIELD_ID]] : u[uOff[C_FIELD_ID]]; 222 const PetscScalar C01 = split ? a[aOff[C_FIELD_ID] + 1] : u[uOff[C_FIELD_ID] + 1]; 223 const PetscScalar C02 = split ? a[aOff[C_FIELD_ID] + 2] : u[uOff[C_FIELD_ID] + 2]; 224 const PetscScalar C11 = split ? a[aOff[C_FIELD_ID] + 3] : u[uOff[C_FIELD_ID] + 3]; 225 const PetscScalar C12 = split ? a[aOff[C_FIELD_ID] + 4] : u[uOff[C_FIELD_ID] + 4]; 226 const PetscScalar C22 = split ? a[aOff[C_FIELD_ID] + 5] : u[uOff[C_FIELD_ID] + 5]; 227 const PetscScalar norm = NORM2C_3d(C00, C01, C02, C11, C12, C22) + eps; 228 const PetscScalar nexp = (gamma - 2.0) / 2.0; 229 const PetscScalar fnorm = PetscPowScalar(norm, nexp); 230 const PetscScalar dfnorm = nexp * PetscPowScalar(norm, nexp - 1.0); 231 const PetscScalar dC[] = {2 * C00, 4 * C01, 4 * C02, 2 * C11, 4 * C12, 2 * C22}; 232 const PetscScalar eqss[] = {0.5, 1., 1., 0.5, 1., 0.5}; 233 234 for (PetscInt k = 0; k < 6; k++) { 235 if (!split) { 236 for (PetscInt j = 0; j < 6; j++) J[k * 6 + j] = eqss[k] * (alpha * dfnorm * dC[j] * u[uOff[C_FIELD_ID] + k]); 237 } 238 J[k * 6 + k] += eqss[k] * (alpha * fnorm + u_tShift); 239 } 240 } 241 242 /* Jacobian for C against C basis functions and gradients of P basis functions */ 243 static void JC_0_c0p1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 244 { 245 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 246 const PetscScalar eqss[] = {0.5, 1., 0.5}; 247 248 J[0] = -2 * gradp[0] * eqss[0]; 249 J[1] = 0.0; 250 251 J[2] = -gradp[1] * eqss[1]; 252 J[3] = -gradp[0] * eqss[1]; 253 254 J[4] = 0.0; 255 J[5] = -2 * gradp[1] * eqss[2]; 256 } 257 258 static void JC_0_c0p1_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 259 { 260 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 261 const PetscScalar eqss[] = {0.5, 1., 1., 0.5, 1., 0.5}; 262 263 J[0] = -2 * gradp[0] * eqss[0]; 264 J[1] = 0.0; 265 J[2] = 0.0; 266 267 J[3] = -gradp[1] * eqss[1]; 268 J[4] = -gradp[0] * eqss[1]; 269 J[5] = 0.0; 270 271 J[6] = -gradp[2] * eqss[2]; 272 J[7] = 0.0; 273 J[8] = -gradp[0] * eqss[2]; 274 275 J[9] = 0.0; 276 J[10] = -2 * gradp[1] * eqss[3]; 277 J[11] = 0.0; 278 279 J[12] = 0.0; 280 J[13] = -gradp[2] * eqss[4]; 281 J[14] = -gradp[1] * eqss[4]; 282 283 J[15] = 0.0; 284 J[16] = 0.0; 285 J[17] = -2 * gradp[2] * eqss[5]; 286 } 287 288 /* residual for C when tested against gradients of basis functions */ 289 static void C_1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 290 { 291 const PetscReal D = PetscRealPart(constants[D_ID]); 292 for (PetscInt k = 0; k < 3; k++) 293 for (PetscInt d = 0; d < 2; d++) f1[k * 2 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 2 + d]; 294 } 295 296 static void C_1_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 297 { 298 const PetscReal D = PetscRealPart(constants[D_ID]); 299 for (PetscInt k = 0; k < 6; k++) 300 for (PetscInt d = 0; d < 3; d++) f1[k * 3 + d] = PetscSqr(D) * u_x[uOff_x[C_FIELD_ID] + k * 3 + d]; 301 } 302 303 /* Jacobian for C against gradients of C basis functions */ 304 static void JC_1_c1c1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 305 { 306 const PetscReal D = PetscRealPart(constants[D_ID]); 307 for (PetscInt k = 0; k < 3; k++) 308 for (PetscInt d = 0; d < 2; d++) J[k * (3 + 1) * 2 * 2 + d * 2 + d] = PetscSqr(D); 309 } 310 311 static void JC_1_c1c1_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 312 { 313 const PetscReal D = PetscRealPart(constants[D_ID]); 314 for (PetscInt k = 0; k < 6; k++) 315 for (PetscInt d = 0; d < 3; d++) J[k * (6 + 1) * 3 * 3 + d * 3 + d] = PetscSqr(D); 316 } 317 318 /* residual for P when tested against basis functions. 319 The source term always comes from the auxiliary data because it must be zero mean (algebraically) */ 320 static void P_0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 321 { 322 PetscScalar S = a[aOff[NUM_FIELDS]]; 323 324 f0[0] = S; 325 } 326 327 /* residual for P when tested against basis functions for the initial condition problem 328 here we don't impose symmetry, and we thus flip the sign of the source function */ 329 static void P_0_aux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 330 { 331 PetscScalar S = a[aOff[NUM_FIELDS]]; 332 333 f0[0] = -S; 334 } 335 336 /* residual for P when tested against gradients of basis functions */ 337 static void P_1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 338 { 339 const PetscReal r = PetscRealPart(constants[R_ID]); 340 const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r; 341 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 342 const PetscScalar C10 = C01; 343 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r; 344 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 345 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0); 346 const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0; 347 348 f1[0] = -((C00 + s) * gradp[0] + C01 * gradp[1]); 349 f1[1] = -(C10 * gradp[0] + (C11 + s) * gradp[1]); 350 } 351 352 static void P_1_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 353 { 354 const PetscReal r = PetscRealPart(constants[R_ID]); 355 const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r; 356 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 357 const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2]; 358 const PetscScalar C10 = C01; 359 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3] + r; 360 const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4]; 361 const PetscScalar C20 = C02; 362 const PetscScalar C21 = C12; 363 const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5] + r; 364 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 365 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0); 366 const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0; 367 368 f1[0] = -((C00 + s) * gradp[0] + C01 * gradp[1] + C02 * gradp[2]); 369 f1[1] = -(C10 * gradp[0] + (C11 + s) * gradp[1] + C12 * gradp[2]); 370 f1[2] = -(C20 * gradp[0] + C21 * gradp[1] + (C22 + s) * gradp[2]); 371 } 372 373 /* Same as above except that the conductivity values come from the auxiliary vec */ 374 static void P_1_aux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 375 { 376 const PetscReal r = PetscRealPart(constants[R_ID]); 377 const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r; 378 const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1]; 379 const PetscScalar C10 = C01; 380 const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2] + r; 381 const PetscScalar *gradp = u_x + uOff_x[Nf > 1 ? P_FIELD_ID : 0]; 382 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0); 383 const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0; 384 385 f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1]; 386 f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1]; 387 } 388 389 static void P_1_aux_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 390 { 391 const PetscReal r = PetscRealPart(constants[R_ID]); 392 const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r; 393 const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1]; 394 const PetscScalar C02 = a[aOff[C_FIELD_ID] + 2]; 395 const PetscScalar C10 = C01; 396 const PetscScalar C11 = a[aOff[C_FIELD_ID] + 3] + r; 397 const PetscScalar C12 = a[aOff[C_FIELD_ID] + 4]; 398 const PetscScalar C20 = C02; 399 const PetscScalar C21 = C12; 400 const PetscScalar C22 = a[aOff[C_FIELD_ID] + 5] + r; 401 const PetscScalar *gradp = u_x + uOff_x[Nf > 1 ? P_FIELD_ID : 0]; 402 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 1.0); 403 const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0; 404 405 f1[0] = (C00 + s) * gradp[0] + C01 * gradp[1] + C02 * gradp[2]; 406 f1[1] = C10 * gradp[0] + (C11 + s) * gradp[1] + C12 * gradp[2]; 407 f1[2] = C20 * gradp[0] + C21 * gradp[1] + (C22 + s) * gradp[2]; 408 } 409 410 /* Jacobian for P against gradients of P basis functions */ 411 static void JP_1_p1p1(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 412 { 413 const PetscReal r = PetscRealPart(constants[R_ID]); 414 const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r; 415 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 416 const PetscScalar C10 = C01; 417 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r; 418 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0); 419 const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0; 420 421 J[0] = -(C00 + s); 422 J[1] = -C01; 423 J[2] = -C10; 424 J[3] = -(C11 + s); 425 } 426 427 static void JP_1_p1p1_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 428 { 429 const PetscReal r = PetscRealPart(constants[R_ID]); 430 const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r; 431 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 432 const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2]; 433 const PetscScalar C10 = C01; 434 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3] + r; 435 const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4]; 436 const PetscScalar C20 = C02; 437 const PetscScalar C21 = C12; 438 const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5] + r; 439 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0); 440 const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0; 441 442 J[0] = -(C00 + s); 443 J[1] = -C01; 444 J[2] = -C02; 445 J[3] = -C10; 446 J[4] = -(C11 + s); 447 J[5] = -C12; 448 J[6] = -C20; 449 J[7] = -C21; 450 J[8] = -(C22 + s); 451 } 452 453 static void JP_1_p1p1_aux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 454 { 455 const PetscReal r = PetscRealPart(constants[R_ID]); 456 const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r; 457 const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1]; 458 const PetscScalar C10 = C01; 459 const PetscScalar C11 = a[aOff[C_FIELD_ID] + 2] + r; 460 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0); 461 const PetscScalar s = fix_c ? FIX_C(C00, C01, C11) : 0.0; 462 463 J[0] = C00 + s; 464 J[1] = C01; 465 J[2] = C10; 466 J[3] = C11 + s; 467 } 468 469 static void JP_1_p1p1_aux_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 470 { 471 const PetscReal r = PetscRealPart(constants[R_ID]); 472 const PetscScalar C00 = a[aOff[C_FIELD_ID]] + r; 473 const PetscScalar C01 = a[aOff[C_FIELD_ID] + 1]; 474 const PetscScalar C02 = a[aOff[C_FIELD_ID] + 2]; 475 const PetscScalar C10 = C01; 476 const PetscScalar C11 = a[aOff[C_FIELD_ID] + 3] + r; 477 const PetscScalar C12 = a[aOff[C_FIELD_ID] + 4]; 478 const PetscScalar C20 = C02; 479 const PetscScalar C21 = C12; 480 const PetscScalar C22 = a[aOff[C_FIELD_ID] + 5] + r; 481 const PetscBool fix_c = (PetscBool)(PetscRealPart(constants[FIXC_ID]) > 0.0); 482 const PetscScalar s = fix_c ? FIX_C_3d(C00, C01, C02, C11, C12, C22) : 0.0; 483 484 J[0] = C00 + s; 485 J[1] = C01; 486 J[2] = C02; 487 J[3] = C10; 488 J[4] = C11 + s; 489 J[5] = C12; 490 J[6] = C20; 491 J[7] = C21; 492 J[8] = C22 + s; 493 } 494 495 /* Jacobian for P against gradients of P basis functions and C basis functions */ 496 static void JP_1_p1c0(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 497 { 498 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 499 500 J[0] = -gradp[0]; 501 J[1] = 0; 502 503 J[2] = -gradp[1]; 504 J[3] = -gradp[0]; 505 506 J[4] = 0; 507 J[5] = -gradp[1]; 508 } 509 510 static void JP_1_p1c0_3d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar J[]) 511 { 512 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 513 514 J[0] = -gradp[0]; 515 J[1] = 0; 516 J[2] = 0; 517 518 J[3] = -gradp[1]; 519 J[4] = -gradp[0]; 520 J[5] = 0; 521 522 J[6] = -gradp[2]; 523 J[7] = 0; 524 J[8] = -gradp[0]; 525 526 J[9] = 0; 527 J[10] = -gradp[1]; 528 J[11] = 0; 529 530 J[12] = 0; 531 J[13] = -gradp[2]; 532 J[14] = -gradp[1]; 533 534 J[15] = 0; 535 J[16] = 0; 536 J[17] = -gradp[2]; 537 } 538 539 /* a collection of gaussian, Dirac-like, source term S(x) = scale * cos(2*pi*(frequency*t + phase)) * exp(-w*||x - x0||^2) */ 540 typedef struct { 541 PetscInt n; 542 PetscReal *x0; 543 PetscReal *w; 544 PetscReal *k; 545 PetscReal *p; 546 PetscReal *r; 547 } MultiSourceCtx; 548 549 typedef struct { 550 PetscReal x0[3]; 551 PetscReal w; 552 PetscReal k; 553 PetscReal p; 554 PetscReal r; 555 } SingleSourceCtx; 556 557 static PetscErrorCode gaussian(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx) 558 { 559 SingleSourceCtx *s = (SingleSourceCtx *)ctx; 560 const PetscReal *x0 = s->x0; 561 const PetscReal w = s->w; 562 const PetscReal k = s->k; /* frequency */ 563 const PetscReal p = s->p; /* phase */ 564 const PetscReal r = s->r; /* scale */ 565 PetscReal n = 0; 566 567 for (PetscInt d = 0; d < dim; ++d) n += (x[d] - x0[d]) * (x[d] - x0[d]); 568 u[0] = r * PetscCosReal(2 * PETSC_PI * (k * time + p)) * PetscExpReal(-w * n); 569 return PETSC_SUCCESS; 570 } 571 572 static PetscErrorCode source(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx) 573 { 574 MultiSourceCtx *s = (MultiSourceCtx *)ctx; 575 576 u[0] = 0.0; 577 for (PetscInt i = 0; i < s->n; i++) { 578 SingleSourceCtx sctx; 579 PetscScalar ut[1]; 580 581 sctx.x0[0] = s->x0[dim * i]; 582 sctx.x0[1] = s->x0[dim * i + 1]; 583 sctx.x0[2] = dim > 2 ? s->x0[dim * i + 2] : 0.0; 584 sctx.w = s->w[i]; 585 sctx.k = s->k[i]; 586 sctx.p = s->p[i]; 587 sctx.r = s->r[i]; 588 589 PetscCall(gaussian(dim, time, x, Nf, ut, &sctx)); 590 591 u[0] += ut[0]; 592 } 593 return PETSC_SUCCESS; 594 } 595 596 /* functionals to be integrated: average -> \int_\Omega u dx */ 597 static void average(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]) 598 { 599 const PetscInt fid = (PetscInt)PetscRealPart(constants[numConstants]); 600 obj[0] = u[uOff[fid]]; 601 } 602 603 /* functionals to be integrated: volume -> \int_\Omega dx */ 604 static void volume(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]) 605 { 606 obj[0] = 1; 607 } 608 609 /* functionals to be integrated: energy -> D^2/2 * ||\nabla C||^2 + c^2\nabla p * (r + C) * \nabla p + \alpha/ \gamma * ||C||^\gamma */ 610 static void energy(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]) 611 { 612 const PetscReal D = PetscRealPart(constants[D_ID]); 613 const PetscReal r = PetscRealPart(constants[R_ID]); 614 const PetscReal alpha = PetscRealPart(constants[ALPHA_ID]); 615 const PetscReal gamma = PetscRealPart(constants[GAMMA_ID]); 616 const PetscReal eps = PetscRealPart(constants[EPS_ID]); 617 618 if (dim == 2) { 619 const PetscScalar C00 = u[uOff[C_FIELD_ID]]; 620 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 621 const PetscScalar C10 = C01; 622 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2]; 623 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 624 const PetscScalar *gradC00 = u_x + uOff_x[C_FIELD_ID]; 625 const PetscScalar *gradC01 = u_x + uOff_x[C_FIELD_ID] + 2; 626 const PetscScalar *gradC11 = u_x + uOff_x[C_FIELD_ID] + 4; 627 const PetscScalar normC = NORM2C(C00, C01, C11) + eps; 628 const PetscScalar normgradC = NORM2C(gradC00[0], gradC01[0], gradC11[0]) + NORM2C(gradC00[1], gradC01[1], gradC11[1]); 629 const PetscScalar nexp = gamma / 2.0; 630 631 const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC; 632 const PetscScalar t1 = gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1]); 633 const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp); 634 635 obj[0] = t0 + t1 + t2; 636 } else { 637 const PetscScalar C00 = u[uOff[C_FIELD_ID]]; 638 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 639 const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2]; 640 const PetscScalar C10 = C01; 641 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3]; 642 const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4]; 643 const PetscScalar C20 = C02; 644 const PetscScalar C21 = C12; 645 const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5]; 646 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 647 const PetscScalar *gradC00 = u_x + uOff_x[C_FIELD_ID]; 648 const PetscScalar *gradC01 = u_x + uOff_x[C_FIELD_ID] + 3; 649 const PetscScalar *gradC02 = u_x + uOff_x[C_FIELD_ID] + 6; 650 const PetscScalar *gradC11 = u_x + uOff_x[C_FIELD_ID] + 9; 651 const PetscScalar *gradC12 = u_x + uOff_x[C_FIELD_ID] + 12; 652 const PetscScalar *gradC22 = u_x + uOff_x[C_FIELD_ID] + 15; 653 const PetscScalar normC = NORM2C_3d(C00, C01, C02, C11, C12, C22) + eps; 654 const PetscScalar normgradC = NORM2C_3d(gradC00[0], gradC01[0], gradC02[0], gradC11[0], gradC12[0], gradC22[0]) + NORM2C_3d(gradC00[1], gradC01[1], gradC02[1], gradC11[1], gradC12[1], gradC22[1]) + NORM2C_3d(gradC00[2], gradC01[2], gradC02[2], gradC11[2], gradC12[2], gradC22[2]); 655 const PetscScalar nexp = gamma / 2.0; 656 657 const PetscScalar t0 = PetscSqr(D) / 2.0 * normgradC; 658 const PetscScalar t1 = gradp[0] * ((C00 + r) * gradp[0] + C01 * gradp[1] + C02 * gradp[2]) + gradp[1] * (C10 * gradp[0] + (C11 + r) * gradp[1] + C12 * gradp[2]) + gradp[2] * (C20 * gradp[0] + C21 * gradp[1] + (C22 + r) * gradp[2]); 659 const PetscScalar t2 = alpha / gamma * PetscPowScalar(normC, nexp); 660 661 obj[0] = t0 + t1 + t2; 662 } 663 } 664 665 /* functionals to be integrated: ellipticity_fail_private -> see below */ 666 static void ellipticity_fail_private(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[], PetscBool add_reg) 667 { 668 #if !PetscDefined(USE_COMPLEX) 669 PetscReal eigs[3]; 670 PetscScalar C00, C01, C02 = 0.0, C11, C12 = 0.0, C22 = 0.0; 671 const PetscReal r = add_reg ? PetscRealPart(constants[R_ID]) : 0.0; 672 673 if (dim == 2) { 674 C00 = u[uOff[C_FIELD_ID]] + r; 675 C01 = u[uOff[C_FIELD_ID] + 1]; 676 C11 = u[uOff[C_FIELD_ID] + 2] + r; 677 } else { 678 C00 = u[uOff[C_FIELD_ID]] + r; 679 C01 = u[uOff[C_FIELD_ID] + 1]; 680 C02 = u[uOff[C_FIELD_ID] + 2]; 681 C11 = u[uOff[C_FIELD_ID] + 3] + r; 682 C12 = u[uOff[C_FIELD_ID] + 4]; 683 C22 = u[uOff[C_FIELD_ID] + 5] + r; 684 } 685 Eigenvalues_Sym3x3(C00, C01, C02, C11, C12, C22, eigs); 686 if (eigs[0] < 0 || eigs[1] < 0 || eigs[2] < 0) obj[0] = -PetscMin(eigs[0], PetscMin(eigs[1], eigs[2])); 687 else obj[0] = 0.0; 688 #else 689 obj[0] = 0.0; 690 #endif 691 } 692 693 /* functionals to be integrated: ellipticity_fail -> 0 means C is elliptic at quadrature point, otherwise it returns 1 */ 694 static void ellipticity_fail(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]) 695 { 696 ellipticity_fail_private(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, obj, PETSC_FALSE); 697 if (PetscAbsScalar(obj[0]) > 0.0) obj[0] = 1.0; 698 } 699 700 /* functionals to be integrated: jacobian_fail -> 0 means C + r is elliptic at quadrature point, otherwise it returns the absolute value of the most negative eigenvalue */ 701 static void jacobian_fail(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]) 702 { 703 ellipticity_fail_private(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, obj, PETSC_TRUE); 704 } 705 706 /* initial conditions for C: eq. 16 */ 707 static PetscErrorCode initial_conditions_C_0(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) 708 { 709 if (dim == 2) { 710 u[0] = 1; 711 u[1] = 0; 712 u[2] = 1; 713 } else { 714 u[0] = 1; 715 u[1] = 0; 716 u[2] = 0; 717 u[3] = 1; 718 u[4] = 0; 719 u[5] = 1; 720 } 721 return PETSC_SUCCESS; 722 } 723 724 /* initial conditions for C: eq. 17 */ 725 static PetscErrorCode initial_conditions_C_1(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) 726 { 727 const PetscReal x = xx[0]; 728 const PetscReal y = xx[1]; 729 730 if (dim == 3) return PETSC_ERR_SUP; 731 u[0] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y)); 732 u[1] = 0; 733 u[2] = (2 - PetscAbsReal(x + y)) * PetscExpReal(-10 * PetscAbsReal(x - y)); 734 return PETSC_SUCCESS; 735 } 736 737 /* initial conditions for C: eq. 18 */ 738 static PetscErrorCode initial_conditions_C_2(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) 739 { 740 u[0] = 0; 741 u[1] = 0; 742 u[2] = 0; 743 if (dim == 3) { 744 u[3] = 0; 745 u[4] = 0; 746 u[5] = 0; 747 } 748 return PETSC_SUCCESS; 749 } 750 751 /* random initial conditions for the diagonal of C */ 752 static PetscErrorCode initial_conditions_C_random(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) 753 { 754 PetscScalar vals[3]; 755 PetscRandom r = (PetscRandom)ctx; 756 757 PetscCall(PetscRandomGetValues(r, dim, vals)); 758 if (dim == 2) { 759 u[0] = vals[0]; 760 u[1] = 0; 761 u[2] = vals[1]; 762 } else { 763 u[0] = vals[0]; 764 u[1] = 0; 765 u[2] = 0; 766 u[3] = vals[1]; 767 u[4] = 0; 768 u[5] = vals[2]; 769 } 770 return PETSC_SUCCESS; 771 } 772 773 /* functionals to be sampled: zero */ 774 static void zero(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) 775 { 776 f[0] = 0.0; 777 } 778 779 /* functionals to be sampled: - (C + r) * \grad p */ 780 static void flux(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) 781 { 782 const PetscReal r = PetscRealPart(constants[R_ID]); 783 const PetscScalar *gradp = u_x + uOff_x[P_FIELD_ID]; 784 785 if (dim == 2) { 786 const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r; 787 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 788 const PetscScalar C10 = C01; 789 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2] + r; 790 791 f[0] = -C00 * gradp[0] - C01 * gradp[1]; 792 f[1] = -C10 * gradp[0] - C11 * gradp[1]; 793 } else { 794 const PetscScalar C00 = u[uOff[C_FIELD_ID]] + r; 795 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 796 const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2]; 797 const PetscScalar C10 = C01; 798 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3] + r; 799 const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4]; 800 const PetscScalar C20 = C02; 801 const PetscScalar C21 = C12; 802 const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5] + r; 803 804 f[0] = -C00 * gradp[0] - C01 * gradp[1] - C02 * gradp[2]; 805 f[1] = -C10 * gradp[0] - C11 * gradp[1] - C12 * gradp[2]; 806 f[2] = -C20 * gradp[0] - C21 * gradp[1] - C22 * gradp[2]; 807 } 808 } 809 810 /* functionals to be sampled: ||C|| */ 811 static void normc(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) 812 { 813 if (dim == 2) { 814 const PetscScalar C00 = u[uOff[C_FIELD_ID]]; 815 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 816 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 2]; 817 818 f[0] = PetscSqrtReal(PetscRealPart(NORM2C(C00, C01, C11))); 819 } else { 820 const PetscScalar C00 = u[uOff[C_FIELD_ID]]; 821 const PetscScalar C01 = u[uOff[C_FIELD_ID] + 1]; 822 const PetscScalar C02 = u[uOff[C_FIELD_ID] + 2]; 823 const PetscScalar C11 = u[uOff[C_FIELD_ID] + 3]; 824 const PetscScalar C12 = u[uOff[C_FIELD_ID] + 4]; 825 const PetscScalar C22 = u[uOff[C_FIELD_ID] + 5]; 826 827 f[0] = PetscSqrtReal(PetscRealPart(NORM2C_3d(C00, C01, C02, C11, C12, C22))); 828 } 829 } 830 831 /* functionals to be sampled: eigenvalues of C */ 832 static void eigsc(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) 833 { 834 #if !PetscDefined(USE_COMPLEX) 835 PetscReal eigs[3]; 836 PetscScalar C00, C01, C02 = 0.0, C11, C12 = 0.0, C22 = 0.0; 837 if (dim == 2) { 838 C00 = u[uOff[C_FIELD_ID]]; 839 C01 = u[uOff[C_FIELD_ID] + 1]; 840 C11 = u[uOff[C_FIELD_ID] + 2]; 841 } else { 842 C00 = u[uOff[C_FIELD_ID]]; 843 C01 = u[uOff[C_FIELD_ID] + 1]; 844 C02 = u[uOff[C_FIELD_ID] + 2]; 845 C11 = u[uOff[C_FIELD_ID] + 3]; 846 C12 = u[uOff[C_FIELD_ID] + 4]; 847 C22 = u[uOff[C_FIELD_ID] + 5]; 848 } 849 Eigenvalues_Sym3x3(C00, C01, C02, C11, C12, C22, eigs); 850 PetscCallVoid(PetscSortReal(dim, eigs)); 851 for (PetscInt d = 0; d < dim; d++) f[d] = eigs[d]; 852 #else 853 for (PetscInt d = 0; d < dim; d++) f[d] = 0; 854 #endif 855 } 856 857 /* functions to be sampled: zero function */ 858 static PetscErrorCode zerof(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) 859 { 860 for (PetscInt d = 0; d < Nc; ++d) u[d] = 0.0; 861 return PETSC_SUCCESS; 862 } 863 864 /* functions to be sampled: constant function */ 865 static PetscErrorCode constantf(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) 866 { 867 for (PetscInt d = 0; d < Nc; ++d) u[d] = 1.0; 868 return PETSC_SUCCESS; 869 } 870 871 /* application context: customizable parameters */ 872 typedef struct { 873 PetscInt dim; 874 PetscBool simplex; 875 PetscReal r; 876 PetscReal eps; 877 PetscReal alpha; 878 PetscReal gamma; 879 PetscReal D; 880 PetscReal domain_volume; 881 PetscInt ic_num; 882 PetscBool split; 883 PetscBool lump; 884 PetscBool amr; 885 PetscBool load; 886 char load_filename[PETSC_MAX_PATH_LEN]; 887 PetscBool save; 888 char save_filename[PETSC_MAX_PATH_LEN]; 889 PetscInt save_every; 890 PetscBool test_restart; 891 PetscInt fix_c; 892 MultiSourceCtx *source_ctx; 893 DM view_dm; 894 TSMonitorVTKCtx view_vtk_ctx; 895 PetscViewerAndFormat *view_hdf5_ctx; 896 PetscInt diagnostic_num; 897 PetscReal view_times[64]; 898 PetscInt view_times_n, view_times_k; 899 PetscReal function_domain_error_tol; 900 VecScatter subsct[NUM_FIELDS]; 901 Vec subvec[NUM_FIELDS]; 902 PetscBool monitor_norms; 903 PetscBool exclude_potential_lte; 904 905 /* hack: need some more plumbing in the library */ 906 SNES snes; 907 } AppCtx; 908 909 /* process command line options */ 910 #include <petsc/private/tsimpl.h> /* To access TSMonitorVTKCtx */ 911 static PetscErrorCode ProcessOptions(AppCtx *options) 912 { 913 char vtkmonfilename[PETSC_MAX_PATH_LEN]; 914 char hdf5monfilename[PETSC_MAX_PATH_LEN]; 915 PetscInt tmp; 916 PetscBool flg; 917 918 PetscFunctionBeginUser; 919 options->dim = 2; 920 options->r = 1.e-1; 921 options->eps = 1.e-3; 922 options->alpha = 0.75; 923 options->gamma = 0.75; 924 options->D = 1.e-2; 925 options->ic_num = 0; 926 options->split = PETSC_FALSE; 927 options->lump = PETSC_FALSE; 928 options->amr = PETSC_FALSE; 929 options->load = PETSC_FALSE; 930 options->save = PETSC_FALSE; 931 options->save_every = -1; 932 options->test_restart = PETSC_FALSE; 933 options->fix_c = 1; /* 1 means only Jac, 2 means function and Jac, < 0 means raise SNESFunctionDomainError when C+r is not posdef */ 934 options->view_vtk_ctx = NULL; 935 options->view_hdf5_ctx = NULL; 936 options->view_dm = NULL; 937 options->diagnostic_num = 1; 938 options->function_domain_error_tol = -1; 939 options->monitor_norms = PETSC_FALSE; 940 options->exclude_potential_lte = PETSC_FALSE; 941 for (PetscInt i = 0; i < NUM_FIELDS; i++) { 942 options->subsct[i] = NULL; 943 options->subvec[i] = NULL; 944 } 945 for (PetscInt i = 0; i < 64; i++) options->view_times[i] = PETSC_MAX_REAL; 946 947 PetscOptionsBegin(PETSC_COMM_WORLD, "", __FILE__, "DMPLEX"); 948 PetscCall(PetscOptionsInt("-dim", "space dimension", __FILE__, options->dim, &options->dim, NULL)); 949 PetscCall(PetscOptionsReal("-alpha", "alpha", __FILE__, options->alpha, &options->alpha, NULL)); 950 PetscCall(PetscOptionsReal("-gamma", "gamma", __FILE__, options->gamma, &options->gamma, NULL)); 951 PetscCall(PetscOptionsReal("-d", "D", __FILE__, options->D, &options->D, NULL)); 952 PetscCall(PetscOptionsReal("-eps", "eps", __FILE__, options->eps, &options->eps, NULL)); 953 PetscCall(PetscOptionsReal("-r", "r", __FILE__, options->r, &options->r, NULL)); 954 PetscCall(PetscOptionsInt("-ic_num", "ic_num", __FILE__, options->ic_num, &options->ic_num, NULL)); 955 PetscCall(PetscOptionsBool("-split", "Operator splitting", __FILE__, options->split, &options->split, NULL)); 956 PetscCall(PetscOptionsBool("-lump", "use mass lumping", __FILE__, options->lump, &options->lump, NULL)); 957 PetscCall(PetscOptionsInt("-fix_c", "Fix conductivity: shift to always be positive semi-definite or raise domain error", __FILE__, options->fix_c, &options->fix_c, NULL)); 958 PetscCall(PetscOptionsBool("-amr", "use adaptive mesh refinement", __FILE__, options->amr, &options->amr, NULL)); 959 PetscCall(PetscOptionsReal("-domain_error_tol", "domain error tolerance", __FILE__, options->function_domain_error_tol, &options->function_domain_error_tol, NULL)); 960 PetscCall(PetscOptionsBool("-test_restart", "test restarting files", __FILE__, options->test_restart, &options->test_restart, NULL)); 961 if (!options->test_restart) { 962 PetscCall(PetscOptionsString("-load", "filename with data to be loaded for restarting", __FILE__, options->load_filename, options->load_filename, PETSC_MAX_PATH_LEN, &options->load)); 963 PetscCall(PetscOptionsString("-save", "filename with data to be saved for restarting", __FILE__, options->save_filename, options->save_filename, PETSC_MAX_PATH_LEN, &options->save)); 964 if (options->save) PetscCall(PetscOptionsInt("-save_every", "save every n timestep (-1 saves only the last)", __FILE__, options->save_every, &options->save_every, NULL)); 965 } 966 PetscCall(PetscOptionsBool("-exclude_potential_lte", "exclude potential from LTE", __FILE__, options->exclude_potential_lte, &options->exclude_potential_lte, NULL)); 967 options->view_times_k = 0; 968 options->view_times_n = 0; 969 PetscCall(PetscOptionsRealArray("-monitor_times", "Save at specific times", NULL, options->view_times, (tmp = 64, &tmp), &flg)); 970 if (flg) options->view_times_n = tmp; 971 972 PetscCall(PetscOptionsString("-monitor_vtk", "Dump VTK file for diagnostic", "TSMonitorSolutionVTK", NULL, vtkmonfilename, sizeof(vtkmonfilename), &flg)); 973 if (flg) { 974 PetscCall(TSMonitorSolutionVTKCtxCreate(vtkmonfilename, &options->view_vtk_ctx)); 975 PetscCall(PetscOptionsInt("-monitor_vtk_interval", "Save every interval time steps", NULL, options->view_vtk_ctx->interval, &options->view_vtk_ctx->interval, NULL)); 976 } 977 PetscCall(PetscOptionsString("-monitor_hdf5", "Dump HDF5 file for diagnostic", "TSMonitorSolution", NULL, hdf5monfilename, sizeof(hdf5monfilename), &flg)); 978 PetscCall(PetscOptionsInt("-monitor_diagnostic_num", "Number of diagnostics to be computed", __FILE__, options->diagnostic_num, &options->diagnostic_num, NULL)); 979 980 if (flg) { 981 #if defined(PETSC_HAVE_HDF5) 982 PetscViewer viewer; 983 984 PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, hdf5monfilename, FILE_MODE_WRITE, &viewer)); 985 PetscCall(PetscViewerAndFormatCreate(viewer, PETSC_VIEWER_HDF5_VIZ, &options->view_hdf5_ctx)); 986 options->view_hdf5_ctx->view_interval = 1; 987 PetscCall(PetscOptionsInt("-monitor_hdf5_interval", "Save every interval time steps", NULL, options->view_hdf5_ctx->view_interval, &options->view_hdf5_ctx->view_interval, NULL)); 988 PetscCall(PetscViewerDestroy(&viewer)); 989 #else 990 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5"); 991 #endif 992 } 993 PetscCall(PetscOptionsBool("-monitor_norms", "Monitor separate SNES norms", __FILE__, options->monitor_norms, &options->monitor_norms, NULL)); 994 995 /* source options */ 996 PetscCall(PetscNew(&options->source_ctx)); 997 options->source_ctx->n = 1; 998 999 PetscCall(PetscOptionsInt("-source_num", "number of sources", __FILE__, options->source_ctx->n, &options->source_ctx->n, NULL)); 1000 tmp = options->source_ctx->n; 1001 PetscCall(PetscMalloc5(options->dim * tmp, &options->source_ctx->x0, tmp, &options->source_ctx->w, tmp, &options->source_ctx->k, tmp, &options->source_ctx->p, tmp, &options->source_ctx->r)); 1002 for (PetscInt i = 0; i < options->source_ctx->n; i++) { 1003 for (PetscInt d = 0; d < options->dim; d++) options->source_ctx->x0[options->dim * i + d] = 0.25; 1004 options->source_ctx->w[i] = 500; 1005 options->source_ctx->k[i] = 0; 1006 options->source_ctx->p[i] = 0; 1007 options->source_ctx->r[i] = 1; 1008 } 1009 tmp = options->dim * options->source_ctx->n; 1010 PetscCall(PetscOptionsRealArray("-source_x0", "source location", __FILE__, options->source_ctx->x0, &tmp, NULL)); 1011 tmp = options->source_ctx->n; 1012 PetscCall(PetscOptionsRealArray("-source_w", "source factor", __FILE__, options->source_ctx->w, &tmp, NULL)); 1013 tmp = options->source_ctx->n; 1014 PetscCall(PetscOptionsRealArray("-source_k", "source frequency", __FILE__, options->source_ctx->k, &tmp, NULL)); 1015 tmp = options->source_ctx->n; 1016 PetscCall(PetscOptionsRealArray("-source_p", "source phase", __FILE__, options->source_ctx->p, &tmp, NULL)); 1017 tmp = options->source_ctx->n; 1018 PetscCall(PetscOptionsRealArray("-source_r", "source scaling", __FILE__, options->source_ctx->r, &tmp, NULL)); 1019 PetscOptionsEnd(); 1020 PetscFunctionReturn(PETSC_SUCCESS); 1021 } 1022 1023 static PetscErrorCode SaveToFile(DM dm, Vec u, const char *filename) 1024 { 1025 #if defined(PETSC_HAVE_HDF5) 1026 PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC; 1027 PetscViewer viewer; 1028 DM cdm = dm; 1029 PetscInt numlevels = 0; 1030 1031 PetscFunctionBeginUser; 1032 while (cdm) { 1033 numlevels++; 1034 PetscCall(DMGetCoarseDM(cdm, &cdm)); 1035 } 1036 /* Cannot be set programmatically */ 1037 PetscCall(PetscOptionsInsertString(NULL, "-dm_plex_view_hdf5_storage_version 3.0.0")); 1038 PetscCall(PetscViewerHDF5Open(PetscObjectComm((PetscObject)dm), filename, FILE_MODE_WRITE, &viewer)); 1039 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "numlevels", PETSC_INT, &numlevels)); 1040 PetscCall(PetscViewerPushFormat(viewer, format)); 1041 for (PetscInt level = numlevels - 1; level >= 0; level--) { 1042 PetscInt cc, rr; 1043 PetscBool isRegular, isUniform; 1044 const char *dmname; 1045 char groupname[PETSC_MAX_PATH_LEN]; 1046 1047 PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level)); 1048 PetscCall(PetscViewerHDF5PushGroup(viewer, groupname)); 1049 PetscCall(PetscObjectGetName((PetscObject)dm, &dmname)); 1050 PetscCall(DMGetCoarsenLevel(dm, &cc)); 1051 PetscCall(DMGetRefineLevel(dm, &rr)); 1052 PetscCall(DMPlexGetRegularRefinement(dm, &isRegular)); 1053 PetscCall(DMPlexGetRefinementUniform(dm, &isUniform)); 1054 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "meshname", PETSC_STRING, dmname)); 1055 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refinelevel", PETSC_INT, &rr)); 1056 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, &cc)); 1057 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refRegular", PETSC_BOOL, &isRegular)); 1058 PetscCall(PetscViewerHDF5WriteAttribute(viewer, NULL, "refUniform", PETSC_BOOL, &isUniform)); 1059 PetscCall(DMPlexTopologyView(dm, viewer)); 1060 PetscCall(DMPlexLabelsView(dm, viewer)); 1061 PetscCall(DMPlexCoordinatesView(dm, viewer)); 1062 PetscCall(DMPlexSectionView(dm, viewer, NULL)); 1063 if (level == numlevels - 1) { 1064 PetscCall(PetscObjectSetName((PetscObject)u, "solution_")); 1065 PetscCall(DMPlexGlobalVectorView(dm, viewer, NULL, u)); 1066 } 1067 if (level) { 1068 PetscInt cStart, cEnd, ccStart, ccEnd, cpStart; 1069 DMPolytopeType ct; 1070 DMPlexTransform tr; 1071 DM sdm; 1072 PetscScalar *array; 1073 PetscSection section; 1074 Vec map; 1075 IS gis; 1076 const PetscInt *gidx; 1077 1078 PetscCall(DMGetCoarseDM(dm, &cdm)); 1079 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1080 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion)); 1081 PetscCall(PetscSectionSetChart(section, cStart, cEnd)); 1082 for (PetscInt c = cStart; c < cEnd; c++) PetscCall(PetscSectionSetDof(section, c, 1)); 1083 PetscCall(PetscSectionSetUp(section)); 1084 1085 PetscCall(DMClone(dm, &sdm)); 1086 PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm")); 1087 PetscCall(PetscObjectSetName((PetscObject)section, "pdm_section")); 1088 PetscCall(DMSetLocalSection(sdm, section)); 1089 PetscCall(PetscSectionDestroy(§ion)); 1090 1091 PetscCall(DMGetLocalVector(sdm, &map)); 1092 PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map")); 1093 PetscCall(VecGetArray(map, &array)); 1094 PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr)); 1095 PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR)); 1096 PetscCall(DMPlexTransformSetDM(tr, cdm)); 1097 PetscCall(DMPlexTransformSetFromOptions(tr)); 1098 PetscCall(DMPlexTransformSetUp(tr)); 1099 PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd)); 1100 PetscCall(DMPlexGetChart(cdm, &cpStart, NULL)); 1101 PetscCall(DMPlexCreatePointNumbering(cdm, &gis)); 1102 PetscCall(ISGetIndices(gis, &gidx)); 1103 for (PetscInt c = ccStart; c < ccEnd; c++) { 1104 PetscInt *rsize, *rcone, *rornt, Nt; 1105 DMPolytopeType *rct; 1106 PetscInt gnum = gidx[c - cpStart] >= 0 ? gidx[c - cpStart] : -(gidx[c - cpStart] + 1); 1107 1108 PetscCall(DMPlexGetCellType(cdm, c, &ct)); 1109 PetscCall(DMPlexTransformCellTransform(tr, ct, c, NULL, &Nt, &rct, &rsize, &rcone, &rornt)); 1110 for (PetscInt r = 0; r < rsize[Nt - 1]; ++r) { 1111 PetscInt pNew; 1112 1113 PetscCall(DMPlexTransformGetTargetPoint(tr, ct, rct[Nt - 1], c, r, &pNew)); 1114 array[pNew - cStart] = gnum; 1115 } 1116 } 1117 PetscCall(ISRestoreIndices(gis, &gidx)); 1118 PetscCall(ISDestroy(&gis)); 1119 PetscCall(VecRestoreArray(map, &array)); 1120 PetscCall(DMPlexTransformDestroy(&tr)); 1121 PetscCall(DMPlexSectionView(dm, viewer, sdm)); 1122 PetscCall(DMPlexLocalVectorView(dm, viewer, sdm, map)); 1123 PetscCall(DMRestoreLocalVector(sdm, &map)); 1124 PetscCall(DMDestroy(&sdm)); 1125 } 1126 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1127 PetscCall(DMGetCoarseDM(dm, &dm)); 1128 } 1129 PetscCall(PetscViewerPopFormat(viewer)); 1130 PetscCall(PetscViewerDestroy(&viewer)); 1131 PetscFunctionReturn(PETSC_SUCCESS); 1132 #else 1133 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5"); 1134 #endif 1135 } 1136 1137 static PetscErrorCode LoadFromFile(MPI_Comm comm, const char *filename, DM *odm) 1138 { 1139 #if defined(PETSC_HAVE_HDF5) 1140 PetscViewerFormat format = PETSC_VIEWER_HDF5_PETSC; 1141 PetscViewer viewer; 1142 DM dm, cdm = NULL; 1143 PetscSF sfXC = NULL; 1144 PetscInt numlevels = -1; 1145 1146 PetscFunctionBeginUser; 1147 PetscCall(PetscViewerHDF5Open(comm, filename, FILE_MODE_READ, &viewer)); 1148 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "numlevels", PETSC_INT, NULL, &numlevels)); 1149 PetscCall(PetscViewerPushFormat(viewer, format)); 1150 for (PetscInt level = 0; level < numlevels; level++) { 1151 char groupname[PETSC_MAX_PATH_LEN], *dmname; 1152 PetscSF sfXB, sfBC, sfG; 1153 PetscPartitioner part; 1154 PetscInt rr, cc; 1155 PetscBool isRegular, isUniform; 1156 1157 PetscCall(DMCreate(comm, &dm)); 1158 PetscCall(DMSetType(dm, DMPLEX)); 1159 PetscCall(PetscSNPrintf(groupname, sizeof(groupname), "level_%" PetscInt_FMT, level)); 1160 PetscCall(PetscViewerHDF5PushGroup(viewer, groupname)); 1161 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "meshname", PETSC_STRING, NULL, &dmname)); 1162 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refinelevel", PETSC_INT, NULL, &rr)); 1163 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "coarsenlevel", PETSC_INT, NULL, &cc)); 1164 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refRegular", PETSC_BOOL, NULL, &isRegular)); 1165 PetscCall(PetscViewerHDF5ReadAttribute(viewer, NULL, "refUniform", PETSC_BOOL, NULL, &isUniform)); 1166 PetscCall(PetscObjectSetName((PetscObject)dm, dmname)); 1167 PetscCall(DMPlexTopologyLoad(dm, viewer, &sfXB)); 1168 PetscCall(DMPlexLabelsLoad(dm, viewer, sfXB)); 1169 PetscCall(DMPlexCoordinatesLoad(dm, viewer, sfXB)); 1170 PetscCall(DMPlexGetPartitioner(dm, &part)); 1171 if (!level) { /* partition the coarse level only */ 1172 PetscCall(PetscPartitionerSetFromOptions(part)); 1173 } else { /* propagate partitioning information from coarser to finer level */ 1174 DM sdm; 1175 Vec map; 1176 PetscSF sf; 1177 PetscLayout clayout; 1178 PetscScalar *array; 1179 PetscInt *cranks_leaf, *cranks_root, *npoints, *points, *ranks, *starts, *gidxs; 1180 PetscInt nparts, cStart, cEnd, nr, ccStart, ccEnd, cpStart, cpEnd; 1181 PetscMPIInt size, rank; 1182 1183 PetscCall(DMClone(dm, &sdm)); 1184 PetscCall(PetscObjectSetName((PetscObject)sdm, "pdm")); 1185 PetscCall(DMPlexSectionLoad(dm, viewer, sdm, sfXB, NULL, &sf)); 1186 PetscCall(DMGetLocalVector(sdm, &map)); 1187 PetscCall(PetscObjectSetName((PetscObject)map, "pdm_map")); 1188 PetscCall(DMPlexLocalVectorLoad(dm, viewer, sdm, sf, map)); 1189 1190 PetscCallMPI(MPI_Comm_size(comm, &size)); 1191 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1192 nparts = size; 1193 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 1194 PetscCall(DMPlexGetHeightStratum(cdm, 0, &ccStart, &ccEnd)); 1195 PetscCall(DMPlexGetChart(cdm, &cpStart, &cpEnd)); 1196 PetscCall(PetscCalloc1(nparts, &npoints)); 1197 PetscCall(PetscMalloc4(cEnd - cStart, &points, cEnd - cStart, &ranks, nparts + 1, &starts, cEnd - cStart, &gidxs)); 1198 PetscCall(PetscSFGetGraph(sfXC, &nr, NULL, NULL, NULL)); 1199 PetscCall(PetscMalloc2(cpEnd - cpStart, &cranks_leaf, nr, &cranks_root)); 1200 for (PetscInt c = 0; c < cpEnd - cpStart; c++) cranks_leaf[c] = rank; 1201 PetscCall(PetscSFReduceBegin(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE)); 1202 PetscCall(PetscSFReduceEnd(sfXC, MPIU_INT, cranks_leaf, cranks_root, MPI_REPLACE)); 1203 1204 PetscCall(VecGetArray(map, &array)); 1205 for (PetscInt c = 0; c < cEnd - cStart; c++) gidxs[c] = (PetscInt)PetscRealPart(array[c]); 1206 PetscCall(VecRestoreArray(map, &array)); 1207 1208 PetscCall(PetscLayoutCreate(comm, &clayout)); 1209 PetscCall(PetscLayoutSetLocalSize(clayout, nr)); 1210 PetscCall(PetscSFSetGraphLayout(sf, clayout, cEnd - cStart, NULL, PETSC_OWN_POINTER, gidxs)); 1211 PetscCall(PetscLayoutDestroy(&clayout)); 1212 1213 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE)); 1214 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, cranks_root, ranks, MPI_REPLACE)); 1215 PetscCall(PetscSFDestroy(&sf)); 1216 PetscCall(PetscFree2(cranks_leaf, cranks_root)); 1217 for (PetscInt c = 0; c < cEnd - cStart; c++) npoints[ranks[c]]++; 1218 1219 starts[0] = 0; 1220 for (PetscInt c = 0; c < nparts; c++) starts[c + 1] = starts[c] + npoints[c]; 1221 for (PetscInt c = 0; c < cEnd - cStart; c++) points[starts[ranks[c]]++] = c; 1222 PetscCall(PetscPartitionerSetType(part, PETSCPARTITIONERSHELL)); 1223 PetscCall(PetscPartitionerShellSetPartition(part, nparts, npoints, points)); 1224 PetscCall(PetscFree(npoints)); 1225 PetscCall(PetscFree4(points, ranks, starts, gidxs)); 1226 PetscCall(DMRestoreLocalVector(sdm, &map)); 1227 PetscCall(DMDestroy(&sdm)); 1228 } 1229 PetscCall(PetscSFDestroy(&sfXC)); 1230 PetscCall(DMPlexDistribute(dm, 0, &sfBC, odm)); 1231 if (*odm) { 1232 PetscCall(DMDestroy(&dm)); 1233 dm = *odm; 1234 *odm = NULL; 1235 PetscCall(PetscObjectSetName((PetscObject)dm, dmname)); 1236 } 1237 if (sfBC) PetscCall(PetscSFCompose(sfXB, sfBC, &sfXC)); 1238 else { 1239 PetscCall(PetscObjectReference((PetscObject)sfXB)); 1240 sfXC = sfXB; 1241 } 1242 PetscCall(PetscSFDestroy(&sfXB)); 1243 PetscCall(PetscSFDestroy(&sfBC)); 1244 PetscCall(DMSetCoarsenLevel(dm, cc)); 1245 PetscCall(DMSetRefineLevel(dm, rr)); 1246 PetscCall(DMPlexSetRegularRefinement(dm, isRegular)); 1247 PetscCall(DMPlexSetRefinementUniform(dm, isUniform)); 1248 PetscCall(DMPlexSectionLoad(dm, viewer, NULL, sfXC, &sfG, NULL)); 1249 if (level == numlevels - 1) { 1250 Vec u; 1251 1252 PetscCall(DMGetNamedGlobalVector(dm, "solution_", &u)); 1253 PetscCall(PetscObjectSetName((PetscObject)u, "solution_")); 1254 PetscCall(DMPlexGlobalVectorLoad(dm, viewer, NULL, sfG, u)); 1255 PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &u)); 1256 } 1257 PetscCall(PetscFree(dmname)); 1258 PetscCall(PetscSFDestroy(&sfG)); 1259 PetscCall(DMSetCoarseDM(dm, cdm)); 1260 PetscCall(DMDestroy(&cdm)); 1261 PetscCall(PetscViewerHDF5PopGroup(viewer)); 1262 cdm = dm; 1263 } 1264 *odm = dm; 1265 PetscCall(PetscViewerPopFormat(viewer)); 1266 PetscCall(PetscViewerDestroy(&viewer)); 1267 PetscCall(PetscSFDestroy(&sfXC)); 1268 PetscFunctionReturn(PETSC_SUCCESS); 1269 #else 1270 SETERRQ(comm, PETSC_ERR_SUP, "Needs HDF5 support. Please reconfigure using --download-hdf5"); 1271 #endif 1272 } 1273 1274 /* 1275 Setup AuxDM: 1276 - project source function and make it zero-mean 1277 - sample frozen fields for operator splitting 1278 */ 1279 static PetscErrorCode ProjectAuxDM(DM dm, PetscReal time, Vec u, AppCtx *ctx) 1280 { 1281 DM dmAux; 1282 Vec la, a; 1283 void *ctxs[NUM_FIELDS + 1]; 1284 PetscScalar vals[NUM_FIELDS + 1]; 1285 VecScatter sctAux; 1286 PetscDS ds; 1287 PetscErrorCode (*funcs[NUM_FIELDS + 1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1288 1289 PetscFunctionBeginUser; 1290 PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &la)); 1291 if (!la) { 1292 PetscFE field; 1293 PetscInt fields[NUM_FIELDS]; 1294 IS is; 1295 Vec tu, ta; 1296 PetscInt dim; 1297 1298 PetscCall(DMClone(dm, &dmAux)); 1299 PetscCall(DMSetNumFields(dmAux, NUM_FIELDS + 1)); 1300 for (PetscInt i = 0; i < NUM_FIELDS; i++) { 1301 PetscCall(DMGetField(dm, i, NULL, (PetscObject *)&field)); 1302 PetscCall(DMSetField(dmAux, i, NULL, (PetscObject)field)); 1303 fields[i] = i; 1304 } 1305 /* PetscFEDuplicate? */ 1306 PetscCall(DMGetDimension(dm, &dim)); 1307 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, ctx->simplex, "p_", -1, &field)); 1308 PetscCall(PetscObjectSetName((PetscObject)field, "source")); 1309 PetscCall(DMSetField(dmAux, NUM_FIELDS, NULL, (PetscObject)field)); 1310 PetscCall(PetscFEDestroy(&field)); 1311 PetscCall(DMCreateDS(dmAux)); 1312 PetscCall(DMCreateSubDM(dmAux, NUM_FIELDS, fields, &is, NULL)); 1313 PetscCall(DMGetGlobalVector(dm, &tu)); 1314 PetscCall(DMGetGlobalVector(dmAux, &ta)); 1315 PetscCall(VecScatterCreate(tu, NULL, ta, is, &sctAux)); 1316 PetscCall(DMRestoreGlobalVector(dm, &tu)); 1317 PetscCall(DMRestoreGlobalVector(dmAux, &ta)); 1318 PetscCall(PetscObjectCompose((PetscObject)dmAux, "scatterAux", (PetscObject)sctAux)); 1319 PetscCall(VecScatterDestroy(&sctAux)); 1320 PetscCall(ISDestroy(&is)); 1321 PetscCall(DMCreateLocalVector(dmAux, &la)); 1322 PetscCall(PetscObjectSetName((PetscObject)la, "auxiliary_")); 1323 PetscCall(DMSetAuxiliaryVec(dm, NULL, 0, 0, la)); 1324 PetscCall(DMDestroy(&dmAux)); 1325 PetscCall(VecDestroy(&la)); 1326 } 1327 if (time == PETSC_MIN_REAL) PetscFunctionReturn(PETSC_SUCCESS); 1328 1329 PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &la)); 1330 PetscCall(VecGetDM(la, &dmAux)); 1331 PetscCall(DMGetDS(dmAux, &ds)); 1332 PetscCall(DMGetGlobalVector(dmAux, &a)); 1333 funcs[C_FIELD_ID] = zerof; 1334 ctxs[C_FIELD_ID] = NULL; 1335 funcs[P_FIELD_ID] = zerof; 1336 ctxs[P_FIELD_ID] = NULL; 1337 funcs[NUM_FIELDS] = source; 1338 ctxs[NUM_FIELDS] = ctx->source_ctx; 1339 PetscCall(DMProjectFunction(dmAux, time, funcs, ctxs, INSERT_ALL_VALUES, a)); 1340 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero)); 1341 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, zero)); 1342 PetscCall(PetscDSSetObjective(ds, NUM_FIELDS, average)); 1343 PetscCall(DMPlexComputeIntegralFEM(dmAux, a, vals, NULL)); 1344 PetscCall(VecShift(a, -vals[NUM_FIELDS] / ctx->domain_volume)); 1345 if (u) { 1346 PetscCall(PetscObjectQuery((PetscObject)dmAux, "scatterAux", (PetscObject *)&sctAux)); 1347 PetscCheck(sctAux, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing scatterAux"); 1348 PetscCall(VecScatterBegin(sctAux, u, a, INSERT_VALUES, SCATTER_FORWARD)); 1349 PetscCall(VecScatterEnd(sctAux, u, a, INSERT_VALUES, SCATTER_FORWARD)); 1350 } 1351 PetscCall(DMGlobalToLocal(dmAux, a, INSERT_VALUES, la)); 1352 PetscCall(VecViewFromOptions(la, NULL, "-aux_view")); 1353 PetscCall(DMRestoreGlobalVector(dmAux, &a)); 1354 PetscFunctionReturn(PETSC_SUCCESS); 1355 } 1356 1357 /* callback for the creation of the potential null space */ 1358 static PetscErrorCode CreatePotentialNullSpace(DM dm, PetscInt ofield, PetscInt nfield, MatNullSpace *nullSpace) 1359 { 1360 Vec vec; 1361 PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {zerof}; 1362 1363 PetscFunctionBeginUser; 1364 funcs[nfield] = constantf; 1365 PetscCall(DMCreateGlobalVector(dm, &vec)); 1366 PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec)); 1367 PetscCall(VecNormalize(vec, NULL)); 1368 PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace)); 1369 /* break ref cycles */ 1370 PetscCall(VecSetDM(vec, NULL)); 1371 PetscCall(VecDestroy(&vec)); 1372 PetscFunctionReturn(PETSC_SUCCESS); 1373 } 1374 1375 static PetscErrorCode DMGetLumpedMass(DM dm, PetscBool local, Vec *lumped_mass) 1376 { 1377 PetscBool has; 1378 1379 PetscFunctionBeginUser; 1380 if (local) { 1381 PetscCall(DMHasNamedLocalVector(dm, "lumped_mass", &has)); 1382 PetscCall(DMGetNamedLocalVector(dm, "lumped_mass", lumped_mass)); 1383 } else { 1384 PetscCall(DMHasNamedGlobalVector(dm, "lumped_mass", &has)); 1385 PetscCall(DMGetNamedGlobalVector(dm, "lumped_mass", lumped_mass)); 1386 } 1387 if (!has) { 1388 Vec w; 1389 IS is; 1390 1391 PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is)); 1392 PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS"); 1393 if (local) { 1394 Vec w2, wg; 1395 1396 PetscCall(DMCreateMassMatrixLumped(dm, &w, NULL)); 1397 PetscCall(DMGetGlobalVector(dm, &wg)); 1398 PetscCall(DMGetLocalVector(dm, &w2)); 1399 PetscCall(VecSet(w2, 0.0)); 1400 PetscCall(VecSet(wg, 1.0)); 1401 PetscCall(VecISSet(wg, is, 0.0)); 1402 PetscCall(DMGlobalToLocal(dm, wg, INSERT_VALUES, w2)); 1403 PetscCall(VecPointwiseMult(w, w, w2)); 1404 PetscCall(DMRestoreGlobalVector(dm, &wg)); 1405 PetscCall(DMRestoreLocalVector(dm, &w2)); 1406 } else { 1407 PetscCall(DMCreateMassMatrixLumped(dm, NULL, &w)); 1408 PetscCall(VecISSet(w, is, 0.0)); 1409 } 1410 PetscCall(VecCopy(w, *lumped_mass)); 1411 PetscCall(VecDestroy(&w)); 1412 } 1413 PetscFunctionReturn(PETSC_SUCCESS); 1414 } 1415 1416 static PetscErrorCode DMRestoreLumpedMass(DM dm, PetscBool local, Vec *lumped_mass) 1417 { 1418 PetscFunctionBeginUser; 1419 if (local) PetscCall(DMRestoreNamedLocalVector(dm, "lumped_mass", lumped_mass)); 1420 else PetscCall(DMRestoreNamedGlobalVector(dm, "lumped_mass", lumped_mass)); 1421 PetscFunctionReturn(PETSC_SUCCESS); 1422 } 1423 1424 /* callbacks for residual and Jacobian */ 1425 static PetscErrorCode DMPlexTSComputeIFunctionFEM_Private(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 1426 { 1427 Vec work, local_lumped_mass; 1428 AppCtx *ctx = (AppCtx *)user; 1429 1430 PetscFunctionBeginUser; 1431 if (ctx->fix_c < 0) { 1432 PetscReal vals[NUM_FIELDS]; 1433 PetscDS ds; 1434 1435 PetscCall(DMGetDS(dm, &ds)); 1436 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, jacobian_fail)); 1437 PetscCall(DMPlexSNESComputeObjectiveFEM(dm, locX, vals, NULL)); 1438 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy)); 1439 if (vals[C_FIELD_ID]) PetscCall(SNESSetFunctionDomainError(ctx->snes)); 1440 } 1441 if (ctx->lump) { 1442 PetscCall(DMGetLumpedMass(dm, PETSC_TRUE, &local_lumped_mass)); 1443 PetscCall(DMGetLocalVector(dm, &work)); 1444 PetscCall(VecSet(work, 0.0)); 1445 PetscCall(DMPlexTSComputeIFunctionFEM(dm, time, locX, work, locF, user)); 1446 PetscCall(VecPointwiseMult(work, locX_t, local_lumped_mass)); 1447 PetscCall(VecAXPY(locF, 1.0, work)); 1448 PetscCall(DMRestoreLocalVector(dm, &work)); 1449 PetscCall(DMRestoreLumpedMass(dm, PETSC_TRUE, &local_lumped_mass)); 1450 } else { 1451 PetscCall(DMPlexTSComputeIFunctionFEM(dm, time, locX, locX_t, locF, user)); 1452 } 1453 PetscFunctionReturn(PETSC_SUCCESS); 1454 } 1455 1456 static PetscErrorCode DMPlexTSComputeIJacobianFEM_Private(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user) 1457 { 1458 Vec lumped_mass, work; 1459 AppCtx *ctx = (AppCtx *)user; 1460 1461 PetscFunctionBeginUser; 1462 if (ctx->lump) { 1463 PetscCall(DMGetLumpedMass(dm, PETSC_FALSE, &lumped_mass)); 1464 PetscCall(DMPlexTSComputeIJacobianFEM(dm, time, locX, locX_t, 0.0, Jac, JacP, user)); 1465 PetscCall(DMGetGlobalVector(dm, &work)); 1466 PetscCall(VecAXPBY(work, X_tShift, 0.0, lumped_mass)); 1467 PetscCall(MatDiagonalSet(JacP, work, ADD_VALUES)); 1468 PetscCall(DMRestoreGlobalVector(dm, &work)); 1469 PetscCall(DMRestoreLumpedMass(dm, PETSC_FALSE, &lumped_mass)); 1470 } else { 1471 PetscCall(DMPlexTSComputeIJacobianFEM(dm, time, locX, locX_t, X_tShift, Jac, JacP, user)); 1472 } 1473 PetscFunctionReturn(PETSC_SUCCESS); 1474 } 1475 1476 /* customize residuals and Jacobians */ 1477 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) 1478 { 1479 PetscDS ds; 1480 PetscInt cdim, dim; 1481 PetscScalar constants[NUM_CONSTANTS]; 1482 PetscScalar vals[NUM_FIELDS]; 1483 PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID}; 1484 Vec dummy; 1485 IS is; 1486 1487 PetscFunctionBeginUser; 1488 constants[R_ID] = ctx->r; 1489 constants[EPS_ID] = ctx->eps; 1490 constants[ALPHA_ID] = ctx->alpha; 1491 constants[GAMMA_ID] = ctx->gamma; 1492 constants[D_ID] = ctx->D; 1493 constants[FIXC_ID] = ctx->fix_c; 1494 constants[SPLIT_ID] = ctx->split; 1495 1496 PetscCall(DMGetDimension(dm, &dim)); 1497 PetscCall(DMGetCoordinateDim(dm, &cdim)); 1498 PetscCheck(dim == 2 || dim == 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Topological dimension must be 2 or 3"); 1499 PetscCheck(dim == ctx->dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Topological dimension mismatch: expected %" PetscInt_FMT ", found %" PetscInt_FMT, dim, ctx->dim); 1500 PetscCheck(cdim == ctx->dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Geometrical dimension mismatch: expected %" PetscInt_FMT ", found %" PetscInt_FMT, cdim, ctx->dim); 1501 PetscCall(DMGetDS(dm, &ds)); 1502 PetscCall(PetscDSSetConstants(ds, NUM_CONSTANTS, constants)); 1503 PetscCall(PetscDSSetImplicit(ds, C_FIELD_ID, PETSC_TRUE)); 1504 PetscCall(PetscDSSetImplicit(ds, P_FIELD_ID, PETSC_TRUE)); 1505 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy)); 1506 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero)); 1507 if (ctx->dim == 2) { 1508 PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0, C_1)); 1509 PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1)); 1510 PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0, NULL, NULL, ctx->D > 0 ? JC_1_c1c1 : NULL)); 1511 if (!ctx->split) { /* we solve a block diagonal system to mimic operator splitting, jacobians will not be correct */ 1512 PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1, NULL, NULL)); 1513 PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0, NULL)); 1514 } 1515 PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1)); 1516 } else { 1517 PetscCall(PetscDSSetResidual(ds, C_FIELD_ID, C_0_3d, C_1_3d)); 1518 PetscCall(PetscDSSetResidual(ds, P_FIELD_ID, P_0, P_1_3d)); 1519 PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, C_FIELD_ID, JC_0_c0c0_3d, NULL, NULL, ctx->D > 0 ? JC_1_c1c1_3d : NULL)); 1520 if (!ctx->split) { 1521 PetscCall(PetscDSSetJacobian(ds, C_FIELD_ID, P_FIELD_ID, NULL, JC_0_c0p1_3d, NULL, NULL)); 1522 PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, C_FIELD_ID, NULL, NULL, JP_1_p1c0_3d, NULL)); 1523 } 1524 PetscCall(PetscDSSetJacobian(ds, P_FIELD_ID, P_FIELD_ID, NULL, NULL, NULL, JP_1_p1p1_3d)); 1525 } 1526 /* Attach potential nullspace */ 1527 PetscCall(DMSetNullSpaceConstructor(dm, P_FIELD_ID, CreatePotentialNullSpace)); 1528 1529 /* Compute domain volume */ 1530 PetscCall(DMGetGlobalVector(dm, &dummy)); 1531 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, volume)); 1532 PetscCall(DMPlexComputeIntegralFEM(dm, dummy, vals, NULL)); 1533 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero)); 1534 PetscCall(DMRestoreGlobalVector(dm, &dummy)); 1535 ctx->domain_volume = PetscRealPart(vals[P_FIELD_ID]); 1536 1537 /* Index sets for potential and conductivities */ 1538 PetscCall(DMCreateSubDM(dm, 1, fields, &is, NULL)); 1539 PetscCall(PetscObjectCompose((PetscObject)dm, "IS conductivity", (PetscObject)is)); 1540 PetscCall(PetscObjectSetName((PetscObject)is, "C")); 1541 PetscCall(ISViewFromOptions(is, NULL, "-is_conductivity_view")); 1542 PetscCall(ISDestroy(&is)); 1543 PetscCall(DMCreateSubDM(dm, 1, fields + 1, &is, NULL)); 1544 PetscCall(PetscObjectSetName((PetscObject)is, "P")); 1545 PetscCall(PetscObjectCompose((PetscObject)dm, "IS potential", (PetscObject)is)); 1546 PetscCall(ISViewFromOptions(is, NULL, "-is_potential_view")); 1547 PetscCall(ISDestroy(&is)); 1548 1549 /* Create auxiliary DM */ 1550 PetscCall(ProjectAuxDM(dm, PETSC_MIN_REAL, NULL, ctx)); 1551 1552 /* Add callbacks */ 1553 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM_Private, ctx)); 1554 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM_Private, ctx)); 1555 /* DMPlexTSComputeBoundary is not needed because we use natural boundary conditions */ 1556 PetscCall(DMTSSetBoundaryLocal(dm, NULL, NULL)); 1557 1558 /* handle nullspace in residual (move it to TSComputeIFunction_DMLocal?) */ 1559 { 1560 MatNullSpace nullsp; 1561 1562 PetscCall(CreatePotentialNullSpace(dm, P_FIELD_ID, P_FIELD_ID, &nullsp)); 1563 PetscCall(PetscObjectCompose((PetscObject)dm, "__dmtsnullspace", (PetscObject)nullsp)); 1564 PetscCall(MatNullSpaceDestroy(&nullsp)); 1565 } 1566 PetscFunctionReturn(PETSC_SUCCESS); 1567 } 1568 1569 /* setup discrete spaces and residuals */ 1570 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) 1571 { 1572 DM cdm = dm; 1573 PetscFE feC, feP; 1574 PetscInt dim; 1575 MPI_Comm comm = PetscObjectComm((PetscObject)dm); 1576 1577 PetscFunctionBeginUser; 1578 PetscCall(DMGetDimension(dm, &dim)); 1579 1580 /* We model Cij with Cij = Cji -> dim*(dim+1)/2 components */ 1581 PetscCall(PetscFECreateDefault(comm, dim, (dim * (dim + 1)) / 2, ctx->simplex, "c_", -1, &feC)); 1582 PetscCall(PetscObjectSetName((PetscObject)feC, "conductivity")); 1583 PetscCall(PetscFECreateDefault(comm, dim, 1, ctx->simplex, "p_", -1, &feP)); 1584 PetscCall(PetscObjectSetName((PetscObject)feP, "potential")); 1585 PetscCall(PetscFECopyQuadrature(feP, feC)); 1586 PetscCall(PetscFEViewFromOptions(feP, NULL, "-view_fe")); 1587 PetscCall(PetscFEViewFromOptions(feC, NULL, "-view_fe")); 1588 1589 PetscCall(DMSetNumFields(dm, 2)); 1590 PetscCall(DMSetField(dm, C_FIELD_ID, NULL, (PetscObject)feC)); 1591 PetscCall(DMSetField(dm, P_FIELD_ID, NULL, (PetscObject)feP)); 1592 PetscCall(PetscFEDestroy(&feC)); 1593 PetscCall(PetscFEDestroy(&feP)); 1594 PetscCall(DMCreateDS(dm)); 1595 while (cdm) { 1596 PetscCall(DMCopyDisc(dm, cdm)); 1597 PetscCall(SetupProblem(cdm, ctx)); 1598 PetscCall(DMGetCoarseDM(cdm, &cdm)); 1599 } 1600 PetscFunctionReturn(PETSC_SUCCESS); 1601 } 1602 1603 /* Create mesh by command line options */ 1604 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) 1605 { 1606 DM plex; 1607 1608 PetscFunctionBeginUser; 1609 if (ctx->load) { 1610 PetscInt refine = 0; 1611 PetscBool isHierarchy; 1612 DM *dms; 1613 char typeName[256]; 1614 PetscBool flg; 1615 1616 PetscCall(LoadFromFile(comm, ctx->load_filename, dm)); 1617 PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX"); 1618 PetscCall(PetscOptionsFList("-dm_mat_type", "Matrix type used for created matrices", "DMSetMatType", MatList, MATAIJ, typeName, sizeof(typeName), &flg)); 1619 if (flg) PetscCall(DMSetMatType(*dm, typeName)); 1620 PetscCall(PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL, 0)); 1621 PetscCall(PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy, 0)); 1622 PetscOptionsEnd(); 1623 if (refine) { 1624 PetscCall(SetupDiscretization(*dm, ctx)); 1625 PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_TRUE)); 1626 } 1627 PetscCall(PetscCalloc1(refine, &dms)); 1628 if (isHierarchy) PetscCall(DMRefineHierarchy(*dm, refine, dms)); 1629 for (PetscInt r = 0; r < refine; r++) { 1630 Mat M; 1631 DM dmr = dms[r]; 1632 Vec u, ur; 1633 1634 if (!isHierarchy) { 1635 PetscCall(DMRefine(*dm, PetscObjectComm((PetscObject)*dm), &dmr)); 1636 PetscCall(DMSetCoarseDM(dmr, *dm)); 1637 } 1638 PetscCall(DMCreateInterpolation(*dm, dmr, &M, NULL)); 1639 PetscCall(DMGetNamedGlobalVector(*dm, "solution_", &u)); 1640 PetscCall(DMGetNamedGlobalVector(dmr, "solution_", &ur)); 1641 PetscCall(MatInterpolate(M, u, ur)); 1642 PetscCall(DMRestoreNamedGlobalVector(*dm, "solution_", &u)); 1643 PetscCall(DMRestoreNamedGlobalVector(dmr, "solution_", &ur)); 1644 PetscCall(MatDestroy(&M)); 1645 if (!isHierarchy) PetscCall(DMSetCoarseDM(dmr, NULL)); 1646 PetscCall(DMDestroy(dm)); 1647 *dm = dmr; 1648 } 1649 if (refine && !isHierarchy) PetscCall(DMSetRefineLevel(*dm, 0)); 1650 PetscCall(PetscFree(dms)); 1651 } else { 1652 PetscCall(DMCreate(comm, dm)); 1653 PetscCall(DMSetType(*dm, DMPLEX)); 1654 PetscCall(DMSetFromOptions(*dm)); 1655 PetscCall(DMLocalizeCoordinates(*dm)); 1656 { 1657 char convType[256]; 1658 PetscBool flg; 1659 PetscOptionsBegin(comm, "", "Additional mesh options", "DMPLEX"); 1660 PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg)); 1661 PetscOptionsEnd(); 1662 if (flg) { 1663 DM dmConv; 1664 PetscCall(DMConvert(*dm, convType, &dmConv)); 1665 if (dmConv) { 1666 PetscCall(DMDestroy(dm)); 1667 *dm = dmConv; 1668 PetscCall(DMSetFromOptions(*dm)); 1669 PetscCall(DMSetUp(*dm)); 1670 } 1671 } 1672 } 1673 } 1674 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, "ref_")); 1675 PetscCall(DMPlexDistributeSetDefault(*dm, PETSC_FALSE)); 1676 PetscCall(DMSetFromOptions(*dm)); 1677 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*dm, NULL)); 1678 1679 PetscCall(DMConvert(*dm, DMPLEX, &plex)); 1680 PetscCall(DMPlexIsSimplex(plex, &ctx->simplex)); 1681 PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->simplex, 1, MPI_C_BOOL, MPI_LOR, comm)); 1682 PetscCall(DMDestroy(&plex)); 1683 1684 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 1685 PetscFunctionReturn(PETSC_SUCCESS); 1686 } 1687 1688 /* Make potential field zero mean */ 1689 static PetscErrorCode ZeroMeanPotential(DM dm, Vec u, PetscScalar domain_volume) 1690 { 1691 PetscScalar vals[NUM_FIELDS]; 1692 PetscDS ds; 1693 IS is; 1694 1695 PetscFunctionBeginUser; 1696 PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is)); 1697 PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS"); 1698 PetscCall(DMGetDS(dm, &ds)); 1699 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average)); 1700 PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL)); 1701 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero)); 1702 PetscCall(VecISShift(u, is, -vals[P_FIELD_ID] / domain_volume)); 1703 PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL)); 1704 PetscFunctionReturn(PETSC_SUCCESS); 1705 } 1706 1707 /* Compute initial conditions and exclude potential from local truncation error 1708 Since we are solving a DAE, once the initial conditions for the differential 1709 variables are set, we need to compute the corresponding value for the 1710 algebraic variables. We do so by creating a subDM for the potential only 1711 and solve a static problem with SNES */ 1712 static PetscErrorCode SetInitialConditionsAndTolerances(TS ts, PetscInt nv, Vec vecs[], PetscBool valid) 1713 { 1714 DM dm; 1715 Vec u, p, subaux, vatol, vrtol; 1716 PetscReal t, atol, rtol; 1717 PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID}; 1718 IS isp; 1719 DM dmp; 1720 VecScatter sctp = NULL; 1721 PetscDS ds; 1722 SNES snes; 1723 KSP ksp; 1724 PC pc; 1725 AppCtx *ctx; 1726 1727 PetscFunctionBeginUser; 1728 PetscCall(TSGetDM(ts, &dm)); 1729 PetscCall(TSGetApplicationContext(ts, &ctx)); 1730 PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&isp)); 1731 PetscCheck(isp, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS"); 1732 if (valid) { 1733 if (ctx->exclude_potential_lte) { 1734 PetscCall(DMCreateGlobalVector(dm, &vatol)); 1735 PetscCall(DMCreateGlobalVector(dm, &vrtol)); 1736 PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL)); 1737 PetscCall(VecSet(vatol, atol)); 1738 PetscCall(VecISSet(vatol, isp, -1)); 1739 PetscCall(VecSet(vrtol, rtol)); 1740 PetscCall(VecISSet(vrtol, isp, -1)); 1741 PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol)); 1742 PetscCall(VecDestroy(&vatol)); 1743 PetscCall(VecDestroy(&vrtol)); 1744 } 1745 for (PetscInt i = 0; i < nv; i++) PetscCall(ZeroMeanPotential(dm, vecs[i], ctx->domain_volume)); 1746 PetscFunctionReturn(PETSC_SUCCESS); 1747 } 1748 PetscCall(DMCreateSubDM(dm, 1, fields + 1, NULL, &dmp)); 1749 PetscCall(DMSetMatType(dmp, MATAIJ)); 1750 PetscCall(DMGetDS(dmp, &ds)); 1751 if (ctx->dim == 2) { 1752 PetscCall(PetscDSSetResidual(ds, 0, P_0_aux, P_1_aux)); 1753 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux)); 1754 } else { 1755 PetscCall(PetscDSSetResidual(ds, 0, P_0_aux, P_1_aux_3d)); 1756 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, JP_1_p1p1_aux_3d)); 1757 } 1758 PetscCall(DMPlexSetSNESLocalFEM(dmp, PETSC_FALSE, NULL)); 1759 1760 PetscCall(DMCreateGlobalVector(dmp, &p)); 1761 1762 PetscCall(SNESCreate(PetscObjectComm((PetscObject)dmp), &snes)); 1763 PetscCall(SNESSetDM(snes, dmp)); 1764 PetscCall(SNESSetOptionsPrefix(snes, "initial_")); 1765 PetscCall(SNESSetErrorIfNotConverged(snes, PETSC_TRUE)); 1766 PetscCall(SNESGetKSP(snes, &ksp)); 1767 PetscCall(KSPSetType(ksp, KSPFGMRES)); 1768 PetscCall(KSPGetPC(ksp, &pc)); 1769 PetscCall(PCSetType(pc, PCGAMG)); 1770 PetscCall(SNESSetFromOptions(snes)); 1771 PetscCall(SNESSetUp(snes)); 1772 1773 /* Loop over input vectors and compute corresponding potential */ 1774 for (PetscInt i = 0; i < nv; i++) { 1775 u = vecs[i]; 1776 if (!valid) { /* Assumes entries in u are not valid */ 1777 PetscErrorCode (*funcs[NUM_FIELDS])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 1778 void *ctxs[NUM_FIELDS] = {NULL}; 1779 PetscRandom r = NULL; 1780 1781 PetscCall(TSGetTime(ts, &t)); 1782 switch (ctx->ic_num) { 1783 case 0: 1784 funcs[C_FIELD_ID] = initial_conditions_C_0; 1785 break; 1786 case 1: 1787 funcs[C_FIELD_ID] = initial_conditions_C_1; 1788 break; 1789 case 2: 1790 funcs[C_FIELD_ID] = initial_conditions_C_2; 1791 break; 1792 case 3: 1793 funcs[C_FIELD_ID] = initial_conditions_C_random; 1794 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)ts), &r)); 1795 PetscCall(PetscRandomSetOptionsPrefix(r, "ic_")); 1796 PetscCall(PetscRandomSetFromOptions(r)); 1797 ctxs[C_FIELD_ID] = r; 1798 break; 1799 default: 1800 SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Unknown IC"); 1801 } 1802 funcs[P_FIELD_ID] = zerof; 1803 PetscCall(DMProjectFunction(dm, t, funcs, ctxs, INSERT_ALL_VALUES, u)); 1804 PetscCall(ProjectAuxDM(dm, t, u, ctx)); 1805 PetscCall(PetscRandomDestroy(&r)); 1806 } 1807 1808 /* pass conductivity information via auxiliary data */ 1809 PetscCall(DMGetAuxiliaryVec(dm, NULL, 0, 0, &subaux)); 1810 PetscCall(DMSetAuxiliaryVec(dmp, NULL, 0, 0, subaux)); 1811 1812 /* solve the subproblem */ 1813 if (!sctp) PetscCall(VecScatterCreate(u, isp, p, NULL, &sctp)); 1814 PetscCall(VecScatterBegin(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD)); 1815 PetscCall(VecScatterEnd(sctp, u, p, INSERT_VALUES, SCATTER_FORWARD)); 1816 PetscCall(SNESSolve(snes, NULL, p)); 1817 1818 /* scatter from potential only to full space */ 1819 PetscCall(VecScatterBegin(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE)); 1820 PetscCall(VecScatterEnd(sctp, p, u, INSERT_VALUES, SCATTER_REVERSE)); 1821 PetscCall(ZeroMeanPotential(dm, u, ctx->domain_volume)); 1822 } 1823 PetscCall(VecDestroy(&p)); 1824 PetscCall(DMDestroy(&dmp)); 1825 PetscCall(SNESDestroy(&snes)); 1826 PetscCall(VecScatterDestroy(&sctp)); 1827 1828 /* exclude potential from computation of the LTE */ 1829 if (ctx->exclude_potential_lte) { 1830 PetscCall(DMCreateGlobalVector(dm, &vatol)); 1831 PetscCall(DMCreateGlobalVector(dm, &vrtol)); 1832 PetscCall(TSGetTolerances(ts, &atol, NULL, &rtol, NULL)); 1833 PetscCall(VecSet(vatol, atol)); 1834 PetscCall(VecISSet(vatol, isp, -1)); 1835 PetscCall(VecSet(vrtol, rtol)); 1836 PetscCall(VecISSet(vrtol, isp, -1)); 1837 PetscCall(TSSetTolerances(ts, atol, vatol, rtol, vrtol)); 1838 PetscCall(VecDestroy(&vatol)); 1839 PetscCall(VecDestroy(&vrtol)); 1840 } 1841 PetscFunctionReturn(PETSC_SUCCESS); 1842 } 1843 1844 /* mesh adaption context */ 1845 typedef struct { 1846 VecTagger refineTag; 1847 DMLabel adaptLabel; 1848 PetscInt cnt; 1849 } AdaptCtx; 1850 1851 static PetscErrorCode ResizeSetUp(TS ts, PetscInt nstep, PetscReal time, Vec u, PetscBool *resize, void *vctx) 1852 { 1853 AdaptCtx *ctx = (AdaptCtx *)vctx; 1854 Vec ellVecCells, ellVecCellsF; 1855 DM dm, plex; 1856 PetscDS ds; 1857 PetscReal norm; 1858 PetscInt cStart, cEnd; 1859 1860 PetscFunctionBeginUser; 1861 PetscCall(TSGetDM(ts, &dm)); 1862 PetscCall(DMConvert(dm, DMPLEX, &plex)); 1863 PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd)); 1864 PetscCall(DMDestroy(&plex)); 1865 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), NUM_FIELDS * (cEnd - cStart), PETSC_DECIDE, &ellVecCellsF)); 1866 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ts), cEnd - cStart, PETSC_DECIDE, &ellVecCells)); 1867 PetscCall(VecSetBlockSize(ellVecCellsF, NUM_FIELDS)); 1868 PetscCall(DMGetDS(dm, &ds)); 1869 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail)); 1870 PetscCall(DMPlexComputeCellwiseIntegralFEM(dm, u, ellVecCellsF, NULL)); 1871 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy)); 1872 PetscCall(VecStrideGather(ellVecCellsF, C_FIELD_ID, ellVecCells, INSERT_VALUES)); 1873 PetscCall(VecDestroy(&ellVecCellsF)); 1874 PetscCall(VecNorm(ellVecCells, NORM_1, &norm)); 1875 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "STEP %d norm %g\n", (int)nstep, (double)norm)); 1876 if (norm && !ctx->cnt) { 1877 IS refineIS; 1878 1879 *resize = PETSC_TRUE; 1880 if (!ctx->refineTag) { 1881 VecTaggerBox refineBox; 1882 refineBox.min = PETSC_MACHINE_EPSILON; 1883 refineBox.max = PETSC_MAX_REAL; 1884 1885 PetscCall(VecTaggerCreate(PETSC_COMM_SELF, &ctx->refineTag)); 1886 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)ctx->refineTag, "refine_")); 1887 PetscCall(VecTaggerSetType(ctx->refineTag, VECTAGGERABSOLUTE)); 1888 PetscCall(VecTaggerAbsoluteSetBox(ctx->refineTag, &refineBox)); 1889 PetscCall(VecTaggerSetFromOptions(ctx->refineTag)); 1890 PetscCall(VecTaggerSetUp(ctx->refineTag)); 1891 PetscCall(PetscObjectViewFromOptions((PetscObject)ctx->refineTag, NULL, "-tag_view")); 1892 } 1893 PetscCall(DMLabelDestroy(&ctx->adaptLabel)); 1894 PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)ts), "adapt", &ctx->adaptLabel)); 1895 PetscCall(VecTaggerComputeIS(ctx->refineTag, ellVecCells, &refineIS, NULL)); 1896 PetscCall(DMLabelSetStratumIS(ctx->adaptLabel, DM_ADAPT_REFINE, refineIS)); 1897 PetscCall(ISDestroy(&refineIS)); 1898 #if 0 1899 void (*funcs[NUM_FIELDS])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]); 1900 Vec ellVec; 1901 1902 funcs[P_FIELD_ID] = ellipticity_fail; 1903 funcs[C_FIELD_ID] = NULL; 1904 1905 PetscCall(DMGetGlobalVector(dm, &ellVec)); 1906 PetscCall(DMProjectField(dm, 0, u, funcs, INSERT_VALUES, ellVec)); 1907 PetscCall(VecViewFromOptions(ellVec,NULL,"-view_amr_ell")); 1908 PetscCall(DMRestoreGlobalVector(dm, &ellVec)); 1909 #endif 1910 ctx->cnt++; 1911 } else { 1912 ctx->cnt = 0; 1913 } 1914 PetscCall(VecDestroy(&ellVecCells)); 1915 PetscFunctionReturn(PETSC_SUCCESS); 1916 } 1917 1918 static PetscErrorCode ResizeTransfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *vctx) 1919 { 1920 AdaptCtx *actx = (AdaptCtx *)vctx; 1921 AppCtx *ctx; 1922 DM dm, adm; 1923 PetscReal time; 1924 PetscInt fields[NUM_FIELDS] = {C_FIELD_ID, P_FIELD_ID}; 1925 IS is; 1926 1927 PetscFunctionBeginUser; 1928 PetscCall(TSGetDM(ts, &dm)); 1929 PetscCheck(actx->adaptLabel, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptLabel"); 1930 PetscCall(DMAdaptLabel(dm, actx->adaptLabel, &adm)); 1931 PetscCheck(adm, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adapted DM"); 1932 PetscCall(TSGetTime(ts, &time)); 1933 PetscCall(DMLabelDestroy(&actx->adaptLabel)); 1934 for (PetscInt i = 0; i < nv; i++) { 1935 PetscCall(DMCreateGlobalVector(adm, &vecsout[i])); 1936 PetscCall(DMForestTransferVec(dm, vecsin[i], adm, vecsout[i], PETSC_TRUE, time)); 1937 } 1938 PetscCall(DMForestSetAdaptivityForest(adm, NULL)); 1939 PetscCall(DMSetCoarseDM(adm, NULL)); 1940 PetscCall(DMSetLocalSection(adm, NULL)); 1941 PetscCall(TSSetDM(ts, adm)); 1942 PetscCall(TSGetTime(ts, &time)); 1943 PetscCall(TSGetApplicationContext(ts, &ctx)); 1944 PetscCall(DMSetNullSpaceConstructor(adm, P_FIELD_ID, CreatePotentialNullSpace)); 1945 PetscCall(DMCreateSubDM(adm, 1, fields, &is, NULL)); 1946 PetscCall(PetscObjectCompose((PetscObject)adm, "IS conductivity", (PetscObject)is)); 1947 PetscCall(ISDestroy(&is)); 1948 PetscCall(DMCreateSubDM(adm, 1, fields + 1, &is, NULL)); 1949 PetscCall(PetscObjectCompose((PetscObject)adm, "IS potential", (PetscObject)is)); 1950 PetscCall(ISDestroy(&is)); 1951 PetscCall(ProjectAuxDM(adm, time, NULL, ctx)); 1952 { 1953 MatNullSpace nullsp; 1954 1955 PetscCall(CreatePotentialNullSpace(adm, P_FIELD_ID, P_FIELD_ID, &nullsp)); 1956 PetscCall(PetscObjectCompose((PetscObject)adm, "__dmtsnullspace", (PetscObject)nullsp)); 1957 PetscCall(MatNullSpaceDestroy(&nullsp)); 1958 } 1959 PetscCall(SetInitialConditionsAndTolerances(ts, nv, vecsout, PETSC_TRUE)); 1960 PetscCall(DMDestroy(&ctx->view_dm)); 1961 for (PetscInt i = 0; i < NUM_FIELDS; i++) { 1962 PetscCall(VecScatterDestroy(&ctx->subsct[i])); 1963 PetscCall(VecDestroy(&ctx->subvec[i])); 1964 } 1965 PetscFunctionReturn(PETSC_SUCCESS); 1966 } 1967 1968 static PetscErrorCode ComputeDiagnostic(Vec u, AppCtx *ctx, Vec *out) 1969 { 1970 DM dm; 1971 PetscInt dim; 1972 PetscFE feFluxC, feNormC, feEigsC; 1973 1974 void (*funcs[])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]) = {normc, eigsc, flux}; 1975 1976 PetscFunctionBeginUser; 1977 if (!ctx->view_dm) { 1978 PetscFE feP; 1979 PetscInt nf = PetscMax(PetscMin(ctx->diagnostic_num, 3), 1); 1980 1981 PetscCall(VecGetDM(u, &dm)); 1982 PetscCall(DMGetDimension(dm, &dim)); 1983 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, 1, ctx->simplex, "normc_", -1, &feNormC)); 1984 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, ctx->simplex, "eigsc_", -1, &feEigsC)); 1985 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, ctx->simplex, "flux_", -1, &feFluxC)); 1986 PetscCall(DMGetField(dm, P_FIELD_ID, NULL, (PetscObject *)&feP)); 1987 PetscCall(PetscFECopyQuadrature(feP, feNormC)); 1988 PetscCall(PetscFECopyQuadrature(feP, feEigsC)); 1989 PetscCall(PetscFECopyQuadrature(feP, feFluxC)); 1990 PetscCall(PetscObjectSetName((PetscObject)feNormC, "normC")); 1991 PetscCall(PetscObjectSetName((PetscObject)feEigsC, "eigsC")); 1992 PetscCall(PetscObjectSetName((PetscObject)feFluxC, "flux")); 1993 1994 PetscCall(DMClone(dm, &ctx->view_dm)); 1995 PetscCall(DMSetNumFields(ctx->view_dm, nf)); 1996 PetscCall(DMSetField(ctx->view_dm, 0, NULL, (PetscObject)feNormC)); 1997 if (nf > 1) PetscCall(DMSetField(ctx->view_dm, 1, NULL, (PetscObject)feEigsC)); 1998 if (nf > 2) PetscCall(DMSetField(ctx->view_dm, 2, NULL, (PetscObject)feFluxC)); 1999 PetscCall(DMCreateDS(ctx->view_dm)); 2000 PetscCall(PetscFEDestroy(&feFluxC)); 2001 PetscCall(PetscFEDestroy(&feNormC)); 2002 PetscCall(PetscFEDestroy(&feEigsC)); 2003 } 2004 PetscCall(DMCreateGlobalVector(ctx->view_dm, out)); 2005 PetscCall(DMProjectField(ctx->view_dm, 0.0, u, funcs, INSERT_VALUES, *out)); 2006 PetscFunctionReturn(PETSC_SUCCESS); 2007 } 2008 2009 static PetscErrorCode MakeScatterAndVec(Vec X, IS is, Vec *Y, VecScatter *sct) 2010 { 2011 PetscInt n; 2012 2013 PetscFunctionBeginUser; 2014 PetscCall(ISGetLocalSize(is, &n)); 2015 PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)X), n, PETSC_DECIDE, Y)); 2016 PetscCall(VecScatterCreate(X, is, *Y, NULL, sct)); 2017 PetscFunctionReturn(PETSC_SUCCESS); 2018 } 2019 2020 static PetscErrorCode FunctionDomainError(TS ts, PetscReal time, Vec X, PetscBool *accept) 2021 { 2022 AppCtx *ctx; 2023 PetscScalar vals[NUM_FIELDS]; 2024 DM dm; 2025 PetscDS ds; 2026 2027 PetscFunctionBeginUser; 2028 *accept = PETSC_TRUE; 2029 PetscCall(TSGetApplicationContext(ts, &ctx)); 2030 if (ctx->function_domain_error_tol < 0) PetscFunctionReturn(PETSC_SUCCESS); 2031 PetscCall(TSGetDM(ts, &dm)); 2032 PetscCall(DMGetDS(dm, &ds)); 2033 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail)); 2034 PetscCall(DMPlexComputeIntegralFEM(dm, X, vals, NULL)); 2035 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy)); 2036 if (PetscAbsScalar(vals[C_FIELD_ID]) > ctx->function_domain_error_tol) *accept = PETSC_FALSE; 2037 if (!*accept) PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "Domain error value %g > %g\n", (double)PetscAbsScalar(vals[C_FIELD_ID]), (double)ctx->function_domain_error_tol)); 2038 PetscFunctionReturn(PETSC_SUCCESS); 2039 } 2040 2041 /* Monitor relevant functionals */ 2042 static PetscErrorCode Monitor(TS ts, PetscInt stepnum, PetscReal time, Vec u, void *vctx) 2043 { 2044 PetscScalar vals[2 * NUM_FIELDS]; 2045 DM dm; 2046 PetscDS ds; 2047 AppCtx *ctx; 2048 PetscBool need_hdf5, need_vtk; 2049 2050 PetscFunctionBeginUser; 2051 PetscCall(TSGetDM(ts, &dm)); 2052 PetscCall(TSGetApplicationContext(ts, &ctx)); 2053 PetscCall(DMGetDS(dm, &ds)); 2054 2055 /* monitor energy and potential average */ 2056 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, average)); 2057 PetscCall(DMPlexComputeIntegralFEM(dm, u, vals, NULL)); 2058 PetscCall(PetscDSSetObjective(ds, P_FIELD_ID, zero)); 2059 2060 /* monitor ellipticity_fail */ 2061 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, ellipticity_fail)); 2062 PetscCall(DMPlexComputeIntegralFEM(dm, u, vals + NUM_FIELDS, NULL)); 2063 PetscCall(PetscDSSetObjective(ds, C_FIELD_ID, energy)); 2064 vals[C_FIELD_ID] /= ctx->domain_volume; 2065 2066 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%4" PetscInt_FMT " TS: time %g, energy %g, intp %g, ell %g\n", stepnum, (double)time, (double)PetscRealPart(vals[C_FIELD_ID]), (double)PetscRealPart(vals[P_FIELD_ID]), (double)PetscRealPart(vals[NUM_FIELDS + C_FIELD_ID]))); 2067 2068 /* monitor diagnostic */ 2069 need_hdf5 = (PetscBool)(ctx->view_hdf5_ctx && ((ctx->view_hdf5_ctx->view_interval > 0 && !(stepnum % ctx->view_hdf5_ctx->view_interval)) || (ctx->view_hdf5_ctx->view_interval && ts->reason))); 2070 need_vtk = (PetscBool)(ctx->view_vtk_ctx && ((ctx->view_vtk_ctx->interval > 0 && !(stepnum % ctx->view_vtk_ctx->interval)) || (ctx->view_vtk_ctx->interval && ts->reason))); 2071 if (ctx->view_times_k < ctx->view_times_n && time >= ctx->view_times[ctx->view_times_k] && time < ctx->view_times[ctx->view_times_k + 1]) { 2072 if (ctx->view_hdf5_ctx) need_hdf5 = PETSC_TRUE; 2073 if (ctx->view_vtk_ctx) need_vtk = PETSC_TRUE; 2074 ctx->view_times_k++; 2075 } 2076 if (need_vtk || need_hdf5) { 2077 Vec diagnostic; 2078 PetscBool dump_dm = (PetscBool)(!!ctx->view_dm); 2079 2080 PetscCall(ComputeDiagnostic(u, ctx, &diagnostic)); 2081 if (need_vtk) { 2082 PetscCall(PetscObjectSetName((PetscObject)diagnostic, "")); 2083 PetscCall(TSMonitorSolutionVTK(ts, stepnum, time, diagnostic, ctx->view_vtk_ctx)); 2084 } 2085 if (need_hdf5) { 2086 DM vdm; 2087 PetscInt seqnum; 2088 2089 PetscCall(VecGetDM(diagnostic, &vdm)); 2090 if (!dump_dm) { 2091 PetscCall(PetscViewerPushFormat(ctx->view_hdf5_ctx->viewer, ctx->view_hdf5_ctx->format)); 2092 PetscCall(DMView(vdm, ctx->view_hdf5_ctx->viewer)); 2093 PetscCall(PetscViewerPopFormat(ctx->view_hdf5_ctx->viewer)); 2094 } 2095 PetscCall(DMGetOutputSequenceNumber(vdm, &seqnum, NULL)); 2096 PetscCall(DMSetOutputSequenceNumber(vdm, seqnum + 1, time)); 2097 PetscCall(PetscObjectSetName((PetscObject)diagnostic, "diagnostic")); 2098 PetscCall(PetscViewerPushFormat(ctx->view_hdf5_ctx->viewer, ctx->view_hdf5_ctx->format)); 2099 PetscCall(VecView(diagnostic, ctx->view_hdf5_ctx->viewer)); 2100 if (ctx->diagnostic_num > 3 || ctx->diagnostic_num < 0) { 2101 PetscCall(DMSetOutputSequenceNumber(dm, seqnum + 1, time)); 2102 PetscCall(VecView(u, ctx->view_hdf5_ctx->viewer)); 2103 } 2104 PetscCall(PetscViewerPopFormat(ctx->view_hdf5_ctx->viewer)); 2105 } 2106 PetscCall(VecDestroy(&diagnostic)); 2107 } 2108 PetscFunctionReturn(PETSC_SUCCESS); 2109 } 2110 2111 /* Save restart information */ 2112 static PetscErrorCode MonitorSave(TS ts, PetscInt steps, PetscReal time, Vec u, void *vctx) 2113 { 2114 DM dm; 2115 AppCtx *ctx = (AppCtx *)vctx; 2116 PetscInt save_every = ctx->save_every; 2117 TSConvergedReason reason; 2118 2119 PetscFunctionBeginUser; 2120 if (!ctx->save) PetscFunctionReturn(PETSC_SUCCESS); 2121 PetscCall(TSGetDM(ts, &dm)); 2122 PetscCall(TSGetConvergedReason(ts, &reason)); 2123 if ((save_every > 0 && steps % save_every == 0) || (save_every == -1 && reason) || save_every < -1) PetscCall(SaveToFile(dm, u, ctx->save_filename)); 2124 PetscFunctionReturn(PETSC_SUCCESS); 2125 } 2126 2127 /* Resample source if time dependent */ 2128 static PetscErrorCode PreStage(TS ts, PetscReal stagetime) 2129 { 2130 AppCtx *ctx; 2131 PetscBool resample, ismatis; 2132 Mat A, P; 2133 2134 PetscFunctionBeginUser; 2135 PetscCall(TSGetApplicationContext(ts, &ctx)); 2136 /* in case we need to call SNESSetFunctionDomainError */ 2137 PetscCall(TSGetSNES(ts, &ctx->snes)); 2138 2139 resample = ctx->split ? PETSC_TRUE : PETSC_FALSE; 2140 for (PetscInt i = 0; i < ctx->source_ctx->n; i++) { 2141 if (ctx->source_ctx->k[i] != 0.0) { 2142 resample = PETSC_TRUE; 2143 break; 2144 } 2145 } 2146 if (resample) { 2147 DM dm; 2148 Vec u = NULL; 2149 2150 PetscCall(TSGetDM(ts, &dm)); 2151 /* In case of a multistage method, we always use the frozen values at the previous time-step */ 2152 if (ctx->split) PetscCall(TSGetSolution(ts, &u)); 2153 PetscCall(ProjectAuxDM(dm, stagetime, u, ctx)); 2154 } 2155 2156 /* element matrices are sparse, ignore zero entries */ 2157 PetscCall(TSGetIJacobian(ts, &A, &P, NULL, NULL)); 2158 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis)); 2159 if (!ismatis) PetscCall(MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 2160 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATIS, &ismatis)); 2161 if (!ismatis) PetscCall(MatSetOption(P, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE)); 2162 2163 /* Set symmetric flag */ 2164 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 2165 PetscCall(MatSetOption(P, MAT_SYMMETRIC, PETSC_TRUE)); 2166 PetscCall(MatSetOption(A, MAT_SYMMETRY_ETERNAL, PETSC_TRUE)); 2167 PetscCall(MatSetOption(P, MAT_SYMMETRY_ETERNAL, PETSC_TRUE)); 2168 PetscFunctionReturn(PETSC_SUCCESS); 2169 } 2170 2171 /* Make potential zero mean after SNES solve */ 2172 static PetscErrorCode PostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y) 2173 { 2174 DM dm; 2175 Vec u = Y[stageindex]; 2176 SNES snes; 2177 PetscInt nits, lits, stepnum; 2178 AppCtx *ctx; 2179 2180 PetscFunctionBeginUser; 2181 PetscCall(TSGetDM(ts, &dm)); 2182 PetscCall(TSGetApplicationContext(ts, &ctx)); 2183 2184 PetscCall(ZeroMeanPotential(dm, u, ctx->domain_volume)); 2185 2186 if (ctx->test_restart) PetscFunctionReturn(PETSC_SUCCESS); 2187 2188 /* monitor linear and nonlinear iterations */ 2189 PetscCall(TSGetStepNumber(ts, &stepnum)); 2190 PetscCall(TSGetSNES(ts, &snes)); 2191 PetscCall(SNESGetIterationNumber(snes, &nits)); 2192 PetscCall(SNESGetLinearSolveIterations(snes, &lits)); 2193 2194 /* if function evals in TSDIRK are zero in the first stage, it is FSAL */ 2195 if (stageindex == 0) { 2196 PetscBool dirk; 2197 PetscInt nf; 2198 2199 PetscCall(PetscObjectTypeCompare((PetscObject)ts, TSDIRK, &dirk)); 2200 PetscCall(SNESGetNumberFunctionEvals(snes, &nf)); 2201 if (dirk && nf == 0) nits = 0; 2202 } 2203 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), " step %" PetscInt_FMT " stage %" PetscInt_FMT " nonlinear its %" PetscInt_FMT ", linear its %" PetscInt_FMT "\n", stepnum, stageindex, nits, lits)); 2204 PetscFunctionReturn(PETSC_SUCCESS); 2205 } 2206 2207 PetscErrorCode MonitorNorms(SNES snes, PetscInt its, PetscReal f, void *vctx) 2208 { 2209 AppCtx *ctx = (AppCtx *)vctx; 2210 Vec F; 2211 DM dm; 2212 PetscReal subnorm[NUM_FIELDS]; 2213 2214 PetscFunctionBeginUser; 2215 PetscCall(SNESGetDM(snes, &dm)); 2216 PetscCall(SNESGetFunction(snes, &F, NULL, NULL)); 2217 if (!ctx->subsct[C_FIELD_ID]) { 2218 IS is; 2219 2220 PetscCall(PetscObjectQuery((PetscObject)dm, "IS conductivity", (PetscObject *)&is)); 2221 PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing conductivity IS"); 2222 PetscCall(MakeScatterAndVec(F, is, &ctx->subvec[C_FIELD_ID], &ctx->subsct[C_FIELD_ID])); 2223 } 2224 if (!ctx->subsct[P_FIELD_ID]) { 2225 IS is; 2226 2227 PetscCall(PetscObjectQuery((PetscObject)dm, "IS potential", (PetscObject *)&is)); 2228 PetscCheck(is, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Missing potential IS"); 2229 PetscCall(MakeScatterAndVec(F, is, &ctx->subvec[P_FIELD_ID], &ctx->subsct[P_FIELD_ID])); 2230 } 2231 PetscCall(VecScatterBegin(ctx->subsct[C_FIELD_ID], F, ctx->subvec[C_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD)); 2232 PetscCall(VecScatterEnd(ctx->subsct[C_FIELD_ID], F, ctx->subvec[C_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD)); 2233 PetscCall(VecScatterBegin(ctx->subsct[P_FIELD_ID], F, ctx->subvec[P_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD)); 2234 PetscCall(VecScatterEnd(ctx->subsct[P_FIELD_ID], F, ctx->subvec[P_FIELD_ID], INSERT_VALUES, SCATTER_FORWARD)); 2235 PetscCall(VecNorm(ctx->subvec[C_FIELD_ID], NORM_2, &subnorm[C_FIELD_ID])); 2236 PetscCall(VecNorm(ctx->subvec[P_FIELD_ID], NORM_2, &subnorm[P_FIELD_ID])); 2237 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), " %3" PetscInt_FMT " SNES Function norms %14.12e, %14.12e\n", its, (double)subnorm[C_FIELD_ID], (double)subnorm[P_FIELD_ID])); 2238 PetscFunctionReturn(PETSC_SUCCESS); 2239 } 2240 2241 static PetscErrorCode Run(MPI_Comm comm, AppCtx *ctx) 2242 { 2243 DM dm; 2244 TS ts; 2245 Vec u; 2246 AdaptCtx *actx; 2247 PetscBool flg; 2248 2249 PetscFunctionBeginUser; 2250 if (ctx->test_restart) { 2251 PetscViewer subviewer; 2252 PetscMPIInt rank; 2253 2254 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 2255 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer)); 2256 if (ctx->load) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d loading from %s\n", rank, ctx->load_filename)); 2257 if (ctx->save) PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d saving to %s\n", rank, ctx->save_filename)); 2258 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer)); 2259 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 2260 } else { 2261 PetscCall(PetscPrintf(comm, "----------------------------\n")); 2262 PetscCall(PetscPrintf(comm, "Simulation parameters:\n")); 2263 PetscCall(PetscPrintf(comm, " dim : %" PetscInt_FMT "\n", ctx->dim)); 2264 PetscCall(PetscPrintf(comm, " r : %g\n", (double)ctx->r)); 2265 PetscCall(PetscPrintf(comm, " eps : %g\n", (double)ctx->eps)); 2266 PetscCall(PetscPrintf(comm, " alpha: %g\n", (double)ctx->alpha)); 2267 PetscCall(PetscPrintf(comm, " gamma: %g\n", (double)ctx->gamma)); 2268 PetscCall(PetscPrintf(comm, " D : %g\n", (double)ctx->D)); 2269 if (ctx->load) PetscCall(PetscPrintf(comm, " load : %s\n", ctx->load_filename)); 2270 else PetscCall(PetscPrintf(comm, " IC : %" PetscInt_FMT "\n", ctx->ic_num)); 2271 PetscCall(PetscPrintf(comm, " snum : %" PetscInt_FMT "\n", ctx->source_ctx->n)); 2272 for (PetscInt i = 0; i < ctx->source_ctx->n; i++) { 2273 const PetscReal *x0 = ctx->source_ctx->x0 + ctx->dim * i; 2274 const PetscReal w = ctx->source_ctx->w[i]; 2275 const PetscReal k = ctx->source_ctx->k[i]; 2276 const PetscReal p = ctx->source_ctx->p[i]; 2277 const PetscReal r = ctx->source_ctx->r[i]; 2278 2279 if (ctx->dim == 2) { 2280 PetscCall(PetscPrintf(comm, " x0 : (%g,%g)\n", (double)x0[0], (double)x0[1])); 2281 } else { 2282 PetscCall(PetscPrintf(comm, " x0 : (%g,%g,%g)\n", (double)x0[0], (double)x0[1], (double)x0[2])); 2283 } 2284 PetscCall(PetscPrintf(comm, " scals: (%g,%g,%g,%g)\n", (double)w, (double)k, (double)p, (double)r)); 2285 } 2286 PetscCall(PetscPrintf(comm, "----------------------------\n")); 2287 } 2288 2289 if (!ctx->test_restart) PetscCall(PetscLogStagePush(SetupStage)); 2290 PetscCall(CreateMesh(comm, &dm, ctx)); 2291 PetscCall(SetupDiscretization(dm, ctx)); 2292 2293 PetscCall(TSCreate(comm, &ts)); 2294 PetscCall(TSSetApplicationContext(ts, ctx)); 2295 2296 PetscCall(TSSetDM(ts, dm)); 2297 if (ctx->test_restart) { 2298 PetscViewer subviewer; 2299 PetscMPIInt rank; 2300 PetscInt level; 2301 2302 PetscCall(DMGetRefineLevel(dm, &level)); 2303 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 2304 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer)); 2305 PetscCall(PetscViewerASCIIPrintf(subviewer, "rank %d DM refinement level %" PetscInt_FMT "\n", rank, level)); 2306 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, comm, &subviewer)); 2307 PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD)); 2308 } 2309 2310 if (ctx->test_restart) PetscCall(TSSetMaxSteps(ts, 1)); 2311 PetscCall(TSSetMaxTime(ts, 10.0)); 2312 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 2313 if (!ctx->test_restart) PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL)); 2314 PetscCall(TSMonitorSet(ts, MonitorSave, ctx, NULL)); 2315 PetscCall(PetscNew(&actx)); 2316 if (ctx->amr) PetscCall(TSSetResize(ts, PETSC_TRUE, ResizeSetUp, ResizeTransfer, actx)); 2317 PetscCall(TSSetPreStage(ts, PreStage)); 2318 PetscCall(TSSetPostStage(ts, PostStage)); 2319 PetscCall(TSSetMaxSNESFailures(ts, -1)); 2320 PetscCall(TSSetFunctionDomainError(ts, FunctionDomainError)); 2321 PetscCall(TSSetFromOptions(ts)); 2322 if (ctx->monitor_norms) { 2323 SNES snes; 2324 2325 PetscCall(TSGetSNES(ts, &snes)); 2326 PetscCall(SNESMonitorSet(snes, MonitorNorms, ctx, NULL)); 2327 } 2328 2329 PetscCall(DMCreateGlobalVector(dm, &u)); 2330 PetscCall(PetscObjectSetName((PetscObject)u, "solution_")); 2331 PetscCall(DMHasNamedGlobalVector(dm, "solution_", &flg)); 2332 if (flg) { /* load from restart file */ 2333 Vec ru; 2334 2335 PetscCall(DMGetNamedGlobalVector(dm, "solution_", &ru)); 2336 PetscCall(VecCopy(ru, u)); 2337 PetscCall(DMRestoreNamedGlobalVector(dm, "solution_", &ru)); 2338 } 2339 PetscCall(SetInitialConditionsAndTolerances(ts, 1, &u, flg)); 2340 PetscCall(TSSetSolution(ts, u)); 2341 PetscCall(VecDestroy(&u)); 2342 PetscCall(DMDestroy(&dm)); 2343 if (!ctx->test_restart) PetscCall(PetscLogStagePop()); 2344 2345 if (!ctx->test_restart) PetscCall(PetscLogStagePush(SolveStage)); 2346 PetscCall(TSSolve(ts, NULL)); 2347 if (!ctx->test_restart) PetscCall(PetscLogStagePop()); 2348 if (ctx->view_vtk_ctx) PetscCall(TSMonitorSolutionVTKDestroy(&ctx->view_vtk_ctx)); 2349 if (ctx->view_hdf5_ctx) PetscCall(PetscViewerAndFormatDestroy(&ctx->view_hdf5_ctx)); 2350 PetscCall(DMDestroy(&ctx->view_dm)); 2351 for (PetscInt i = 0; i < NUM_FIELDS; i++) { 2352 PetscCall(VecScatterDestroy(&ctx->subsct[i])); 2353 PetscCall(VecDestroy(&ctx->subvec[i])); 2354 } 2355 2356 PetscCall(TSDestroy(&ts)); 2357 PetscCall(VecTaggerDestroy(&actx->refineTag)); 2358 PetscCall(PetscFree(actx)); 2359 PetscFunctionReturn(PETSC_SUCCESS); 2360 } 2361 2362 int main(int argc, char **argv) 2363 { 2364 AppCtx ctx; 2365 2366 PetscFunctionBeginUser; 2367 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 2368 PetscCall(ProcessOptions(&ctx)); 2369 PetscCall(PetscLogStageRegister("Setup", &SetupStage)); 2370 PetscCall(PetscLogStageRegister("Solve", &SolveStage)); 2371 if (ctx.test_restart) { /* Test sequences of save and loads */ 2372 PetscMPIInt rank; 2373 2374 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 2375 2376 /* test saving */ 2377 ctx.load = PETSC_FALSE; 2378 ctx.save = PETSC_TRUE; 2379 /* sequential save */ 2380 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential save\n")); 2381 PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_seq_%d.h5", rank)); 2382 PetscCall(Run(PETSC_COMM_SELF, &ctx)); 2383 /* parallel save */ 2384 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel save\n")); 2385 PetscCall(PetscSNPrintf(ctx.save_filename, sizeof(ctx.save_filename), "test_ex30_par.h5")); 2386 PetscCall(Run(PETSC_COMM_WORLD, &ctx)); 2387 2388 /* test loading */ 2389 ctx.load = PETSC_TRUE; 2390 ctx.save = PETSC_FALSE; 2391 /* sequential and parallel runs from sequential save */ 2392 PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_seq_0.h5")); 2393 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from sequential save\n")); 2394 PetscCall(Run(PETSC_COMM_SELF, &ctx)); 2395 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from sequential save\n")); 2396 PetscCall(Run(PETSC_COMM_WORLD, &ctx)); 2397 /* sequential and parallel runs from parallel save */ 2398 PetscCall(PetscSNPrintf(ctx.load_filename, sizeof(ctx.load_filename), "test_ex30_par.h5")); 2399 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test sequential load from parallel save\n")); 2400 PetscCall(Run(PETSC_COMM_SELF, &ctx)); 2401 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test parallel load from parallel save\n")); 2402 PetscCall(Run(PETSC_COMM_WORLD, &ctx)); 2403 } else { /* Run the simulation */ 2404 PetscCall(Run(PETSC_COMM_WORLD, &ctx)); 2405 } 2406 PetscCall(PetscFree5(ctx.source_ctx->x0, ctx.source_ctx->w, ctx.source_ctx->k, ctx.source_ctx->p, ctx.source_ctx->r)); 2407 PetscCall(PetscFree(ctx.source_ctx)); 2408 PetscCall(PetscFinalize()); 2409 return 0; 2410 } 2411 2412 /*TEST 2413 2414 testset: 2415 args: -dm_plex_box_faces 3,3 -ksp_type preonly -pc_type svd -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 1 -initial_snes_test_jacobian -snes_test_jacobian -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 3 2416 2417 test: 2418 suffix: 0 2419 nsize: {{1 2}} 2420 args: -dm_refine 1 -lump {{0 1}} -exclude_potential_lte 2421 2422 test: 2423 suffix: 0_split 2424 nsize: {{1 2}} 2425 args: -dm_refine 1 -split 2426 2427 test: 2428 suffix: 0_3d 2429 nsize: {{1 2}} 2430 args: -dm_plex_box_faces 3,3,3 -dim 3 -dm_plex_dim 3 -lump {{0 1}} 2431 2432 test: 2433 suffix: 0_dirk 2434 nsize: {{1 2}} 2435 args: -dm_refine 1 -ts_type dirk 2436 2437 test: 2438 suffix: 0_dirk_mg 2439 nsize: {{1 2}} 2440 args: -dm_refine_hierarchy 1 -ts_type dirk -pc_type mg -mg_levels_pc_type bjacobi -mg_levels_sub_pc_factor_levels 2 -mg_levels_sub_pc_factor_mat_ordering_type rcm -mg_levels_sub_pc_factor_reuse_ordering -mg_coarse_pc_type svd -lump {{0 1}} 2441 2442 test: 2443 suffix: 0_dirk_fieldsplit 2444 nsize: {{1 2}} 2445 args: -dm_refine 1 -ts_type dirk -pc_type fieldsplit -pc_fieldsplit_type multiplicative -lump {{0 1}} 2446 2447 test: 2448 requires: p4est 2449 suffix: 0_p4est 2450 nsize: {{1 2}} 2451 args: -dm_refine 1 -dm_plex_convert_type p4est -lump 0 2452 2453 test: 2454 suffix: 0_periodic 2455 nsize: {{1 2}} 2456 args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -lump {{0 1}} 2457 2458 test: 2459 requires: p4est 2460 suffix: 0_p4est_periodic 2461 nsize: {{1 2}} 2462 args: -dm_plex_box_bd periodic,periodic -dm_refine_pre 1 -dm_plex_convert_type p4est -lump 0 2463 2464 test: 2465 requires: p4est 2466 suffix: 0_p4est_mg 2467 nsize: {{1 2}} 2468 args: -dm_forest_minimum_refinement 0 -dm_forest_initial_refinement 2 -dm_plex_convert_type p4est -pc_type mg -mg_coarse_pc_type svd -mg_levels_pc_type svd -lump 0 2469 2470 testset: 2471 requires: hdf5 2472 args: -test_restart -dm_plex_box_faces 3,3 -ksp_type preonly -pc_type mg -mg_levels_pc_type svd -c_petscspace_degree 1 -p_petscspace_degree 1 -petscpartitioner_type simple -test_restart 2473 2474 test: 2475 requires: !single 2476 suffix: restart 2477 nsize: {{1 2}separate output} 2478 args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 0 2479 2480 test: 2481 requires: triangle 2482 suffix: restart_simplex 2483 nsize: {{1 2}separate output} 2484 args: -dm_refine_hierarchy {{0 1}separate output} -dm_plex_simplex 1 2485 2486 test: 2487 requires: !single 2488 suffix: restart_refonly 2489 nsize: {{1 2}separate output} 2490 args: -dm_refine 1 -dm_plex_simplex 0 2491 2492 test: 2493 requires: triangle 2494 suffix: restart_simplex_refonly 2495 nsize: {{1 2}separate output} 2496 args: -dm_refine 1 -dm_plex_simplex 1 2497 2498 test: 2499 suffix: annulus 2500 requires: exodusii 2501 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -ksp_type preonly -pc_type none -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 2 -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -source_num 2 -source_k 1,2 2502 2503 test: 2504 suffix: hdf5_diagnostic 2505 requires: hdf5 exodusii 2506 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -ksp_type preonly -pc_type none -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 2 -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -source_num 2 -source_k 1,2 -monitor_hdf5 diagnostic.h5 -ic_num 3 2507 2508 test: 2509 suffix: vtk_diagnostic 2510 requires: exodusii 2511 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -ksp_type preonly -pc_type none -c_petscspace_degree 1 -p_petscspace_degree 1 -ts_max_steps 2 -initial_snes_type ksponly -snes_type ksponly -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -source_num 2 -source_k 1,2 -monitor_vtk 'diagnostic-%03d.vtu' -ic_num 3 2512 2513 test: 2514 suffix: full_cdisc 2515 args: -dm_plex_box_faces 3,3 -c_petscspace_degree 0 -p_petscspace_degree 1 -ts_max_steps 1 -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 0 -dm_refine 1 -ts_type beuler -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition selfp -fieldsplit_conductivity_pc_type pbjacobi -fieldsplit_potential_mat_schur_complement_ainv_type blockdiag -fieldsplit_potential_ksp_type preonly -fieldsplit_potential_pc_type svd 2516 2517 test: 2518 suffix: full_cdisc_split 2519 args: -dm_plex_box_faces 3,3 -c_petscspace_degree 0 -p_petscspace_degree 1 -ts_max_steps 1 -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 0 -dm_refine 1 -ts_type beuler -pc_type fieldsplit -pc_fieldsplit_type additive -fieldsplit_conductivity_pc_type pbjacobi -fieldsplit_potential_pc_type gamg -split -monitor_norms 2520 2521 test: 2522 suffix: full_cdisc_minres 2523 args: -dm_plex_box_faces 3,3 -c_petscspace_degree 0 -p_petscspace_degree 1 -ts_max_steps 1 -petscpartitioner_type simple -dm_plex_simplex 0 -ts_adapt_type none -ic_num 0 -dm_refine 1 -ts_type beuler -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type diag -pc_fieldsplit_schur_precondition selfp -fieldsplit_conductivity_pc_type pbjacobi -fieldsplit_potential_mat_schur_complement_ainv_type blockdiag -fieldsplit_potential_ksp_type preonly -fieldsplit_potential_pc_type svd -ksp_type minres 2524 2525 TEST*/ 2526