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 PetscInt cStart, cEnd, cEndInterior; 225 DM plex; 226 PetscErrorCode ierr; 227 228 PetscFunctionBegin; 229 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 230 ierr = DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);CHKERRQ(ierr); 231 ierr = DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 232 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 233 ierr = DMPlexComputeJacobian_Internal(plex, cStart, cEnd, time, X_tShift, locX, locX_t, Jac, JacP, user);CHKERRQ(ierr); 234 ierr = DMDestroy(&plex);CHKERRQ(ierr); 235 PetscFunctionReturn(0); 236 } 237 238 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs) 239 { 240 DM dm; 241 SNES snes; 242 Vec sol; 243 PetscBool check; 244 PetscErrorCode ierr; 245 246 PetscFunctionBegin; 247 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr); 248 if (!check) PetscFunctionReturn(0); 249 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 250 ierr = TSSetSolution(ts, sol);CHKERRQ(ierr); 251 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 252 ierr = TSSetUp(ts);CHKERRQ(ierr); 253 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 254 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr); 255 ierr = DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);CHKERRQ(ierr); 256 ierr = VecDestroy(&sol);CHKERRQ(ierr); 257 PetscFunctionReturn(0); 258 } 259