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