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 = DMPlexGetDataFVM(plex, NULL, &cellgeom, &facegeom, NULL);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 } 80 if (facegeom) {PetscValidPointer(facegeom, 2); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject *) facegeom);CHKERRQ(ierr);} 81 if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject *) cellgeom);CHKERRQ(ierr);} 82 if (minRadius) {ierr = DMPlexGetMinRadius(plex, minRadius);CHKERRQ(ierr);} 83 ierr = DMDestroy(&plex);CHKERRQ(ierr); 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "DMPlexTSGetGradientDM" 89 /*@C 90 DMPlexTSGetGradientDM - Return gradient data layout 91 92 Input Parameters: 93 + dm - The DM 94 - fv - The PetscFV 95 96 Output Parameter: 97 . dmGrad - The layout for gradient values 98 99 Level: developer 100 101 .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal() 102 @*/ 103 PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad) 104 { 105 DMTS dmts; 106 PetscObject obj; 107 PetscBool computeGradients; 108 DM plex; 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 113 PetscValidHeaderSpecific(fv,PETSCFV_CLASSID,2); 114 PetscValidPointer(dmGrad,3); 115 ierr = PetscFVGetComputeGradients(fv, &computeGradients);CHKERRQ(ierr); 116 if (!computeGradients) {*dmGrad = NULL; PetscFunctionReturn(0);} 117 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 118 ierr = DMGetDMTS(plex, &dmts);CHKERRQ(ierr); 119 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);CHKERRQ(ierr); 120 if (!obj) { 121 DM dmGrad; 122 Vec faceGeometry, cellGeometry; 123 124 ierr = DMPlexGetDataFVM(plex,fv,&cellGeometry,&faceGeometry,&dmGrad);CHKERRQ(ierr); 125 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);CHKERRQ(ierr); 126 } 127 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);CHKERRQ(ierr); 128 ierr = DMDestroy(&plex);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131 132 #undef __FUNCT__ 133 #define __FUNCT__ "DMPlexTSComputeRHSFunctionFVM" 134 /*@ 135 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user 136 137 Input Parameters: 138 + dm - The mesh 139 . t - The time 140 . locX - Local solution 141 - user - The user context 142 143 Output Parameter: 144 . F - Global output vector 145 146 Level: developer 147 148 .seealso: DMPlexComputeJacobianActionFEM() 149 @*/ 150 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user) 151 { 152 Vec locF; 153 PetscInt cStart, cEnd, cEndInterior; 154 DM plex; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 159 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 160 ierr = DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 161 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 162 ierr = DMGetLocalVector(plex, &locF);CHKERRQ(ierr); 163 ierr = VecZeroEntries(locF);CHKERRQ(ierr); 164 ierr = DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, NULL, time, locF, user);CHKERRQ(ierr); 165 ierr = DMLocalToGlobalBegin(plex, locF, INSERT_VALUES, F);CHKERRQ(ierr); 166 ierr = DMLocalToGlobalEnd(plex, locF, INSERT_VALUES, F);CHKERRQ(ierr); 167 ierr = DMRestoreLocalVector(plex, &locF);CHKERRQ(ierr); 168 ierr = DMDestroy(&plex);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "DMPlexTSComputeBoundary" 174 /*@ 175 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 176 177 Input Parameters: 178 + dm - The mesh 179 . t - The time 180 . locX - Local solution 181 . locX_t - Local solution time derivative, or NULL 182 - user - The user context 183 184 Level: developer 185 186 .seealso: DMPlexComputeJacobianActionFEM() 187 @*/ 188 PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user) 189 { 190 DM plex; 191 Vec faceGeometryFVM = NULL; 192 PetscInt Nf, f; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 ierr = DMTSConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 197 ierr = DMGetNumFields(plex, &Nf);CHKERRQ(ierr); 198 if (!locX_t) { 199 /* This is the RHS part */ 200 for (f = 0; f < Nf; f++) { 201 PetscObject obj; 202 PetscClassId id; 203 204 ierr = DMGetField(plex, f, &obj);CHKERRQ(ierr); 205 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 206 if (id == PETSCFV_CLASSID) { 207 ierr = DMPlexSNESGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);CHKERRQ(ierr); 208 break; 209 } 210 } 211 } 212 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);CHKERRQ(ierr); 213 /* TODO: locX_t */ 214 ierr = DMDestroy(&plex);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 218 #undef __FUNCT__ 219 #define __FUNCT__ "DMPlexTSComputeIFunctionFEM" 220 /*@ 221 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 222 223 Input Parameters: 224 + dm - The mesh 225 . t - The time 226 . locX - Local solution 227 . locX_t - Local solution time derivative, or NULL 228 - user - The user context 229 230 Output Parameter: 231 . locF - Local output vector 232 233 Level: developer 234 235 .seealso: DMPlexComputeJacobianActionFEM() 236 @*/ 237 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 238 { 239 PetscInt cStart, cEnd, cEndInterior; 240 DM plex; 241 PetscErrorCode ierr; 242 243 PetscFunctionBegin; 244 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 245 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 246 ierr = DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 247 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 248 ierr = DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, locX_t, time, locF, user);CHKERRQ(ierr); 249 ierr = DMDestroy(&plex);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "DMPlexTSComputeIJacobianFEM" 255 /*@ 256 DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user 257 258 Input Parameters: 259 + dm - The mesh 260 . t - The time 261 . locX - Local solution 262 . locX_t - Local solution time derivative, or NULL 263 . X_tshift - The multiplicative parameter for dF/du_t 264 - user - The user context 265 266 Output Parameter: 267 . locF - Local output vector 268 269 Level: developer 270 271 .seealso: DMPlexComputeJacobianActionFEM() 272 @*/ 273 PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user) 274 { 275 PetscInt cStart, cEnd, cEndInterior; 276 DM plex; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 281 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 282 ierr = DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 283 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 284 ierr = DMPlexComputeJacobian_Internal(plex, cStart, cEnd, time, X_tShift, locX, locX_t, Jac, JacP, user);CHKERRQ(ierr); 285 ierr = DMDestroy(&plex);CHKERRQ(ierr); 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "DMTSCheckFromOptions" 291 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs) 292 { 293 DM dm; 294 SNES snes; 295 Vec sol; 296 PetscBool check; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr); 301 if (!check) PetscFunctionReturn(0); 302 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 303 ierr = TSSetSolution(ts, sol);CHKERRQ(ierr); 304 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 305 ierr = TSSetUp(ts);CHKERRQ(ierr); 306 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 307 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr); 308 ierr = DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);CHKERRQ(ierr); 309 ierr = VecDestroy(&sol);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312