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__ "DMPlexTSGetGeometryFVM" 9 /*@ 10 DMPlexTSGetGeometryFVM - Return precomputed geometric data 11 12 Input Parameter: 13 . dm - The DM 14 15 Output Parameters: 16 + facegeom - The values precomputed from face geometry 17 . cellgeom - The values precomputed from cell geometry 18 - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell 19 20 Level: developer 21 22 .seealso: DMPlexTSSetRHSFunctionLocal() 23 @*/ 24 PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius) 25 { 26 DMTS dmts; 27 PetscObject obj; 28 PetscBool isPlex; 29 DM plex = NULL; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 34 ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);CHKERRQ(ierr); 35 if (!isPlex) { 36 ierr = DMConvert(dm,DMPLEX,&plex);CHKERRQ(ierr); 37 if (!plex) {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot convert to DMPlex");} 38 ierr = DMCopyDMTS(dm,plex);CHKERRQ(ierr); 39 dm = plex; 40 } 41 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 42 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", &obj);CHKERRQ(ierr); 43 if (!obj) { 44 Vec cellgeom, facegeom; 45 46 ierr = DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);CHKERRQ(ierr); 47 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject) facegeom);CHKERRQ(ierr); 48 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject) cellgeom);CHKERRQ(ierr); 49 ierr = VecDestroy(&facegeom);CHKERRQ(ierr); 50 ierr = VecDestroy(&cellgeom);CHKERRQ(ierr); 51 } 52 if (facegeom) {PetscValidPointer(facegeom, 2); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject *) facegeom);CHKERRQ(ierr);} 53 if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject *) cellgeom);CHKERRQ(ierr);} 54 if (minRadius) {ierr = DMPlexGetMinRadius(dm, minRadius);CHKERRQ(ierr);} 55 if (!isPlex) { 56 ierr = DMDestroy(&plex);CHKERRQ(ierr); 57 } 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "DMPlexTSGetGradientDM" 63 /*@C 64 DMPlexTSGetGradientDM - Return gradient data layout 65 66 Input Parameters: 67 + dm - The DM 68 - fv - The PetscFV 69 70 Output Parameter: 71 . dmGrad - The layout for gradient values 72 73 Level: developer 74 75 .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal() 76 @*/ 77 PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad) 78 { 79 DMTS dmts; 80 PetscObject obj; 81 PetscBool computeGradients; 82 PetscBool isPlex; 83 DM plex = NULL; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 88 PetscValidHeaderSpecific(fv,PETSCFV_CLASSID,2); 89 PetscValidPointer(dmGrad,3); 90 ierr = PetscFVGetComputeGradients(fv, &computeGradients);CHKERRQ(ierr); 91 if (!computeGradients) {*dmGrad = NULL; PetscFunctionReturn(0);} 92 ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);CHKERRQ(ierr); 93 if (!isPlex) { 94 ierr = DMConvert(dm,DMPLEX,&plex);CHKERRQ(ierr); 95 if (!plex) {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot convert to DMPlex");} 96 ierr = DMCopyDMTS(dm,plex);CHKERRQ(ierr); 97 dm = plex; 98 } 99 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 100 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);CHKERRQ(ierr); 101 if (!obj) { 102 DM dmGrad; 103 Vec faceGeometry, cellGeometry; 104 105 ierr = DMPlexTSGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 106 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);CHKERRQ(ierr); 107 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);CHKERRQ(ierr); 108 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 109 } 110 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);CHKERRQ(ierr); 111 if (!isPlex) { 112 ierr = DMDestroy(&plex);CHKERRQ(ierr); 113 } 114 PetscFunctionReturn(0); 115 } 116 117 #undef __FUNCT__ 118 #define __FUNCT__ "DMPlexTSComputeRHSFunctionFVM" 119 /*@ 120 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user 121 122 Input Parameters: 123 + dm - The mesh 124 . t - The time 125 . locX - Local solution 126 - user - The user context 127 128 Output Parameter: 129 . F - Global output vector 130 131 Level: developer 132 133 .seealso: DMPlexComputeJacobianActionFEM() 134 @*/ 135 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user) 136 { 137 Vec locF; 138 PetscInt cStart, cEnd, cEndInterior; 139 PetscBool isPlex; 140 DM plex = NULL; 141 PetscErrorCode ierr; 142 143 PetscFunctionBegin; 144 ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);CHKERRQ(ierr); 145 if (!isPlex) { 146 ierr = DMConvert(dm,DMPLEX,&plex);CHKERRQ(ierr); 147 if (!plex) {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot convert to DMPlex");} 148 ierr = DMCopyDMTS(dm,plex);CHKERRQ(ierr); 149 dm = plex; 150 } 151 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 152 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 153 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 154 ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 155 ierr = VecZeroEntries(locF);CHKERRQ(ierr); 156 ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, NULL, locF, user);CHKERRQ(ierr); 157 ierr = DMLocalToGlobalBegin(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr); 158 ierr = DMLocalToGlobalEnd(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr); 159 ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 160 if (!isPlex) { 161 ierr = DMDestroy(&plex);CHKERRQ(ierr); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "DMPlexTSComputeIFunctionFEM" 168 /*@ 169 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 170 171 Input Parameters: 172 + dm - The mesh 173 . t - The time 174 . locX - Local solution 175 . locX_t - Local solution time derivative, or NULL 176 - user - The user context 177 178 Output Parameter: 179 . locF - Local output vector 180 181 Level: developer 182 183 .seealso: DMPlexComputeJacobianActionFEM() 184 @*/ 185 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 186 { 187 PetscInt cStart, cEnd, cEndInterior; 188 PetscBool isPlex; 189 DM plex = NULL; 190 PetscErrorCode ierr; 191 192 PetscFunctionBegin; 193 ierr = PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);CHKERRQ(ierr); 194 if (!isPlex) { 195 ierr = DMConvert(dm,DMPLEX,&plex);CHKERRQ(ierr); 196 if (!plex) {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot convert to DMPlex");} 197 ierr = DMCopyDMTS(dm,plex);CHKERRQ(ierr); 198 dm = plex; 199 } 200 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 201 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 202 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 203 ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, locX_t, locF, user);CHKERRQ(ierr); 204 if (!isPlex) { 205 ierr = DMDestroy(&plex);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "DMTSCheckFromOptions" 212 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs) 213 { 214 DM dm; 215 SNES snes; 216 Vec sol; 217 PetscBool check; 218 PetscErrorCode ierr; 219 220 PetscFunctionBegin; 221 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr); 222 if (!check) PetscFunctionReturn(0); 223 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 224 ierr = TSSetSolution(ts, sol);CHKERRQ(ierr); 225 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 226 ierr = TSSetUp(ts);CHKERRQ(ierr); 227 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 228 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr); 229 ierr = DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);CHKERRQ(ierr); 230 ierr = VecDestroy(&sol);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233