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 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user 43 44 Input Parameters: 45 + dm - The mesh 46 . t - The time 47 . locX - Local solution 48 - user - The user context 49 50 Output Parameter: 51 . F - Global output vector 52 53 Level: developer 54 55 .seealso: DMPlexComputeJacobianActionFEM() 56 @*/ 57 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user) 58 { 59 Vec locF; 60 IS cellIS; 61 DM plex; 62 PetscInt depth; 63 PetscErrorCode ierr; 64 65 PetscFunctionBegin; 66 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 67 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 68 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 69 if (!cellIS) { 70 ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr); 71 } 72 ierr = DMGetLocalVector(plex, &locF);CHKERRQ(ierr); 73 ierr = VecZeroEntries(locF);CHKERRQ(ierr); 74 ierr = DMPlexComputeResidual_Internal(plex, cellIS, time, locX, NULL, time, locF, user);CHKERRQ(ierr); 75 ierr = DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);CHKERRQ(ierr); 76 ierr = DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);CHKERRQ(ierr); 77 ierr = DMRestoreLocalVector(plex, &locF);CHKERRQ(ierr); 78 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 79 ierr = DMDestroy(&plex);CHKERRQ(ierr); 80 PetscFunctionReturn(0); 81 } 82 83 /*@ 84 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 85 86 Input Parameters: 87 + dm - The mesh 88 . t - The time 89 . locX - Local solution 90 . locX_t - Local solution time derivative, or NULL 91 - user - The user context 92 93 Level: developer 94 95 .seealso: DMPlexComputeJacobianActionFEM() 96 @*/ 97 PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user) 98 { 99 DM plex; 100 Vec faceGeometryFVM = NULL; 101 PetscInt Nf, f; 102 PetscErrorCode ierr; 103 104 PetscFunctionBegin; 105 ierr = DMTSConvertPlex(dm, &plex, PETSC_TRUE);CHKERRQ(ierr); 106 ierr = DMGetNumFields(plex, &Nf);CHKERRQ(ierr); 107 if (!locX_t) { 108 /* This is the RHS part */ 109 for (f = 0; f < Nf; f++) { 110 PetscObject obj; 111 PetscClassId id; 112 113 ierr = DMGetField(plex, f, NULL, &obj);CHKERRQ(ierr); 114 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 115 if (id == PETSCFV_CLASSID) { 116 ierr = DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);CHKERRQ(ierr); 117 break; 118 } 119 } 120 } 121 ierr = DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);CHKERRQ(ierr); 122 /* TODO: locX_t */ 123 ierr = DMDestroy(&plex);CHKERRQ(ierr); 124 PetscFunctionReturn(0); 125 } 126 127 /*@ 128 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X 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 Output Parameter: 138 . locF - Local output vector 139 140 Level: developer 141 142 .seealso: DMPlexComputeJacobianActionFEM() 143 @*/ 144 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 145 { 146 DM plex; 147 IS cellIS; 148 PetscInt depth; 149 PetscErrorCode ierr; 150 151 PetscFunctionBegin; 152 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 153 ierr = DMPlexGetDepth(plex, &depth);CHKERRQ(ierr); 154 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 155 if (!cellIS) { 156 ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr); 157 } 158 ierr = DMPlexComputeResidual_Internal(plex, cellIS, time, locX, locX_t, time, locF, user);CHKERRQ(ierr); 159 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 160 ierr = DMDestroy(&plex);CHKERRQ(ierr); 161 PetscFunctionReturn(0); 162 } 163 164 /*@ 165 DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user 166 167 Input Parameters: 168 + dm - The mesh 169 . t - The time 170 . locX - Local solution 171 . locX_t - Local solution time derivative, or NULL 172 . X_tshift - The multiplicative parameter for dF/du_t 173 - user - The user context 174 175 Output Parameter: 176 . locF - Local output vector 177 178 Level: developer 179 180 .seealso: DMPlexComputeJacobianActionFEM() 181 @*/ 182 PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user) 183 { 184 DM plex; 185 PetscDS prob; 186 PetscBool hasJac, hasPrec; 187 IS cellIS; 188 PetscInt depth; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 ierr = DMTSConvertPlex(dm,&plex,PETSC_TRUE);CHKERRQ(ierr); 193 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 194 ierr = PetscDSHasJacobian(prob, &hasJac);CHKERRQ(ierr); 195 ierr = PetscDSHasJacobianPreconditioner(prob, &hasPrec);CHKERRQ(ierr); 196 if (hasJac && hasPrec) {ierr = MatZeroEntries(Jac);CHKERRQ(ierr);} 197 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 198 ierr = DMPlexGetDepth(plex,&depth);CHKERRQ(ierr); 199 ierr = DMGetStratumIS(plex, "dim", depth, &cellIS);CHKERRQ(ierr); 200 if (!cellIS) {ierr = DMGetStratumIS(plex, "depth", depth, &cellIS);CHKERRQ(ierr);} 201 ierr = DMPlexComputeJacobian_Internal(plex, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);CHKERRQ(ierr); 202 ierr = ISDestroy(&cellIS);CHKERRQ(ierr); 203 ierr = DMDestroy(&plex);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 /*@C 208 DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information 209 210 Input Parameters: 211 + ts - the TS object 212 - u - representative TS vector 213 214 Note: The user must call PetscDSSetExactSolution() beforehand 215 216 Level: developer 217 @*/ 218 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u) 219 { 220 DM dm; 221 SNES snes; 222 Vec sol; 223 PetscBool check; 224 PetscErrorCode ierr; 225 226 PetscFunctionBegin; 227 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr); 228 if (!check) PetscFunctionReturn(0); 229 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 230 ierr = VecCopy(u, sol);CHKERRQ(ierr); 231 ierr = TSSetSolution(ts, u);CHKERRQ(ierr); 232 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 233 ierr = TSSetUp(ts);CHKERRQ(ierr); 234 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 235 ierr = SNESSetSolution(snes, u);CHKERRQ(ierr); 236 ierr = DMSNESCheck_Internal(snes, dm, sol);CHKERRQ(ierr); 237 ierr = VecDestroy(&sol);CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240