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