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