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