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__ "DMPlexTSGetGeometryFVM" 9 /*@ 10 DMPlexTSGetGeometryFVM - Return precomputed geometric data 11 12 Input Parameter: 13 . dm - The DM 14 15 Output Parameters: 16 + facegeom - The values precomputed from face geometry 17 . cellgeom - The values precomputed from cell geometry 18 - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell 19 20 Level: developer 21 22 .seealso: DMPlexTSSetRHSFunctionLocal() 23 @*/ 24 PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius) 25 { 26 DMTS dmts; 27 PetscObject obj; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 32 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 33 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", &obj);CHKERRQ(ierr); 34 if (!obj) { 35 Vec cellgeom, facegeom; 36 37 ierr = DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);CHKERRQ(ierr); 38 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject) facegeom);CHKERRQ(ierr); 39 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject) cellgeom);CHKERRQ(ierr); 40 ierr = VecDestroy(&facegeom);CHKERRQ(ierr); 41 ierr = VecDestroy(&cellgeom);CHKERRQ(ierr); 42 } 43 if (facegeom) {PetscValidPointer(facegeom, 2); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject *) facegeom);CHKERRQ(ierr);} 44 if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject *) cellgeom);CHKERRQ(ierr);} 45 if (minRadius) {ierr = DMPlexGetMinRadius(dm, minRadius);CHKERRQ(ierr);} 46 PetscFunctionReturn(0); 47 } 48 49 #undef __FUNCT__ 50 #define __FUNCT__ "DMPlexTSGetGradientDM" 51 /*@C 52 DMPlexTSGetGradientDM - Return gradient data layout 53 54 Input Parameters: 55 + dm - The DM 56 - fv - The PetscFV 57 58 Output Parameter: 59 . dmGrad - The layout for gradient values 60 61 Level: developer 62 63 .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal() 64 @*/ 65 PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad) 66 { 67 DMTS dmts; 68 PetscObject obj; 69 PetscBool computeGradients; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 74 PetscValidHeaderSpecific(fv,PETSCFV_CLASSID,2); 75 PetscValidPointer(dmGrad,3); 76 ierr = PetscFVGetComputeGradients(fv, &computeGradients);CHKERRQ(ierr); 77 if (!computeGradients) {*dmGrad = NULL; PetscFunctionReturn(0);} 78 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 79 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);CHKERRQ(ierr); 80 if (!obj) { 81 DM dmGrad; 82 Vec faceGeometry, cellGeometry; 83 84 ierr = DMPlexTSGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 85 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);CHKERRQ(ierr); 86 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);CHKERRQ(ierr); 87 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 88 } 89 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);CHKERRQ(ierr); 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "DMPlexTSComputeRHSFunctionFVM" 95 /*@ 96 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user 97 98 Input Parameters: 99 + dm - The mesh 100 . t - The time 101 . locX - Local solution 102 - user - The user context 103 104 Output Parameter: 105 . F - Global output vector 106 107 Level: developer 108 109 .seealso: DMPlexComputeJacobianActionFEM() 110 @*/ 111 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user) 112 { 113 Vec locF; 114 PetscInt cStart, cEnd, cEndInterior; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 119 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 120 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 121 ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 122 ierr = VecZeroEntries(locF);CHKERRQ(ierr); 123 ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, NULL, locF, user);CHKERRQ(ierr); 124 ierr = DMLocalToGlobalBegin(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr); 125 ierr = DMLocalToGlobalEnd(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr); 126 ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "DMPlexTSComputeIFunctionFEM" 132 /*@ 133 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 134 135 Input Parameters: 136 + dm - The mesh 137 . t - The time 138 . locX - Local solution 139 . locX_t - Local solution time derivative, or NULL 140 - user - The user context 141 142 Output Parameter: 143 . locF - Local output vector 144 145 Level: developer 146 147 .seealso: DMPlexComputeJacobianActionFEM() 148 @*/ 149 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 150 { 151 PetscInt cStart, cEnd, cEndInterior; 152 PetscErrorCode ierr; 153 154 PetscFunctionBegin; 155 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 156 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 157 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 158 ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, locX_t, locF, user);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "DMTSCheckFromOptions" 164 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs) 165 { 166 DM dm; 167 SNES snes; 168 Vec sol; 169 PetscBool check; 170 PetscErrorCode ierr; 171 172 PetscFunctionBegin; 173 ierr = PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);CHKERRQ(ierr); 174 if (!check) PetscFunctionReturn(0); 175 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 176 ierr = TSSetSolution(ts, sol);CHKERRQ(ierr); 177 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 178 ierr = TSSetUp(ts);CHKERRQ(ierr); 179 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 180 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr); 181 ierr = DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);CHKERRQ(ierr); 182 ierr = VecDestroy(&sol);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185