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