1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petsc/private/snesimpl.h> 4 #include <petscds.h> 5 #include <petscfv.h> 6 7 static PetscErrorCode DMTSConvertPlex(DM dm, DM *plex, PetscBool copy) 8 { 9 PetscBool isPlex; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);CHKERRQ(ierr); 14 if (isPlex) { 15 *plex = dm; 16 ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr); 17 } else { 18 ierr = PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);CHKERRQ(ierr); 19 if (!*plex) { 20 ierr = DMConvert(dm,DMPLEX,plex);CHKERRQ(ierr); 21 ierr = PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);CHKERRQ(ierr); 22 if (copy) { 23 PetscInt i; 24 PetscObject obj; 25 const char *comps[3] = {"A","dmAux","dmCh"}; 26 27 ierr = DMCopyDMTS(dm, *plex);CHKERRQ(ierr); 28 ierr = DMCopyDMSNES(dm, *plex);CHKERRQ(ierr); 29 for (i = 0; i < 3; i++) { 30 ierr = PetscObjectQuery((PetscObject) dm, comps[i], &obj);CHKERRQ(ierr); 31 ierr = PetscObjectCompose((PetscObject) *plex, comps[i], obj);CHKERRQ(ierr); 32 } 33 } 34 } else { 35 ierr = PetscObjectReference((PetscObject) *plex);CHKERRQ(ierr); 36 } 37 } 38 PetscFunctionReturn(0); 39 } 40 41 42 /*@ 43 DMPlexTSGetGeometryFVM - Return precomputed geometric data 44 45 Input Parameter: 46 . dm - The DM 47 48 Output Parameters: 49 + facegeom - The values precomputed from face geometry 50 . cellgeom - The values precomputed from cell geometry 51 - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell 52 53 Level: developer 54 55 .seealso: DMPlexTSSetRHSFunctionLocal() 56 @*/ 57 PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius) 58 { 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = DMPlexSNESGetGeometryFVM(dm,facegeom,cellgeom,minRadius);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 /*@C 67 DMPlexTSGetGradientDM - Return gradient data layout 68 69 Input Parameters: 70 + dm - The DM 71 - fv - The PetscFV 72 73 Output Parameter: 74 . dmGrad - The layout for gradient values 75 76 Level: developer 77 78 .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal() 79 @*/ 80 PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 ierr = DMPlexSNESGetGradientDM(dm,fv,dmGrad);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 /*@ 90 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user 91 92 Input Parameters: 93 + dm - The mesh 94 . t - The time 95 . locX - Local solution 96 - user - The user context 97 98 Output Parameter: 99 . F - Global output vector 100 101 Level: developer 102 103 .seealso: DMPlexComputeJacobianActionFEM() 104 @*/ 105 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user) 106 { 107 Vec locF; 108 IS cellIS; 109 DM plex; 110 PetscInt depth; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 115 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 116 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 117 if (!cellIS) { 118 ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr); 119 } 120 ierr = DMGetLocalVector(plex, &locF);CHKERRQ(ierr); 121 ierr = VecZeroEntries(locF);CHKERRQ(ierr); 122 ierr = DMPlexComputeResidual_Internal(plex, cellIS, time, locX, NULL, time, locF, user);CHKERRQ(ierr); 123 ierr = DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);CHKERRQ(ierr); 124 ierr = DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);CHKERRQ(ierr); 125 ierr = DMRestoreLocalVector(plex, &locF);CHKERRQ(ierr); 126 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 127 ierr = DMDestroy(&plex);CHKERRQ(ierr); 128 PetscFunctionReturn(0); 129 } 130 131 /*@ 132 DMPlexTSComputeBoundary - Insert the essential boundary values for the local input X and/or its time derivative X_t using pointwise functions specified by the user 133 134 Input Parameters: 135 + dm - The mesh 136 . t - The time 137 . locX - Local solution 138 . locX_t - Local solution time derivative, or NULL 139 - user - The user context 140 141 Level: developer 142 143 .seealso: DMPlexComputeJacobianActionFEM() 144 @*/ 145 PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user) 146 { 147 DM plex; 148 Vec faceGeometryFVM = NULL; 149 PetscInt Nf, f; 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 ierr = DMTSConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 154 ierr = DMGetNumFields(plex, &Nf);CHKERRQ(ierr); 155 if (!locX_t) { 156 /* This is the RHS part */ 157 for (f = 0; f < Nf; f++) { 158 PetscObject obj; 159 PetscClassId id; 160 161 ierr = DMGetField(plex, f, NULL, &obj);CHKERRQ(ierr); 162 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 163 if (id == PETSCFV_CLASSID) { 164 ierr = DMPlexSNESGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);CHKERRQ(ierr); 165 break; 166 } 167 } 168 } 169 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);CHKERRQ(ierr); 170 /* TODO: locX_t */ 171 ierr = DMDestroy(&plex);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 /*@ 176 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 177 178 Input Parameters: 179 + dm - The mesh 180 . t - The time 181 . locX - Local solution 182 . locX_t - Local solution time derivative, or NULL 183 - user - The user context 184 185 Output Parameter: 186 . locF - Local output vector 187 188 Level: developer 189 190 .seealso: DMPlexComputeJacobianActionFEM() 191 @*/ 192 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 193 { 194 DM plex; 195 IS cellIS; 196 PetscInt depth; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 201 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 202 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 203 if (!cellIS) { 204 ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr); 205 } 206 ierr = DMPlexComputeResidual_Internal(plex, cellIS, time, locX, locX_t, time, locF, user);CHKERRQ(ierr); 207 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 208 ierr = DMDestroy(&plex);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 /*@ 213 DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user 214 215 Input Parameters: 216 + dm - The mesh 217 . t - The time 218 . locX - Local solution 219 . locX_t - Local solution time derivative, or NULL 220 . X_tshift - The multiplicative parameter for dF/du_t 221 - user - The user context 222 223 Output Parameter: 224 . locF - Local output vector 225 226 Level: developer 227 228 .seealso: DMPlexComputeJacobianActionFEM() 229 @*/ 230 PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user) 231 { 232 DM plex; 233 PetscDS prob; 234 PetscBool hasJac, hasPrec; 235 IS cellIS; 236 PetscInt depth; 237 PetscErrorCode ierr; 238 239 PetscFunctionBegin; 240 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 241 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 242 ierr = PetscDSHasJacobian(prob, &hasJac);CHKERRQ(ierr); 243 ierr = PetscDSHasJacobianPreconditioner(prob, &hasPrec);CHKERRQ(ierr); 244 if (hasJac && hasPrec) {ierr = MatZeroEntries(Jac);CHKERRQ(ierr);} 245 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 246 ierr = DMPlexGetDepth(plex,&depth);CHKERRQ(ierr); 247 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 248 if (!cellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);} 249 ierr = DMPlexComputeJacobian_Internal(plex, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);CHKERRQ(ierr); 250 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 251 ierr = DMDestroy(&plex);CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 /*@C 256 DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information 257 258 Input Parameters: 259 + ts - the TS object 260 . u - representative TS vector 261 . exactFuncs - pointwise functions of the exact solution for each field 262 - ctxs - contexts for the functions 263 264 Level: developer 265 @*/ 266 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs) 267 { 268 DM dm; 269 SNES snes; 270 Vec sol; 271 PetscBool check; 272 PetscErrorCode ierr; 273 274 PetscFunctionBegin; 275 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr); 276 if (!check) PetscFunctionReturn(0); 277 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 278 ierr = VecCopy(u, sol);CHKERRQ(ierr); 279 ierr = TSSetSolution(ts, u);CHKERRQ(ierr); 280 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 281 ierr = TSSetUp(ts);CHKERRQ(ierr); 282 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 283 ierr = SNESSetSolution(snes, u);CHKERRQ(ierr); 284 ierr = DMSNESCheck_Internal(snes, dm, sol, exactFuncs, ctxs);CHKERRQ(ierr); 285 ierr = VecDestroy(&sol);CHKERRQ(ierr); 286 PetscFunctionReturn(0); 287 } 288