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