1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscds.h> 4 #include <petscfv.h> 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMPlexTSGetGeometryFVM" 8 /*@ 9 DMPlexTSGetGeometryFVM - Return precomputed geometric data 10 11 Input Parameter: 12 . dm - The DM 13 14 Output Parameters: 15 + facegeom - The values precomputed from face geometry 16 . cellgeom - The values precomputed from cell geometry 17 - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell 18 19 Level: developer 20 21 .seealso: DMPlexTSSetRHSFunctionLocal() 22 @*/ 23 PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius) 24 { 25 DMTS dmts; 26 PetscObject obj; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 31 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 32 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", &obj);CHKERRQ(ierr); 33 if (!obj) { 34 Vec cellgeom, facegeom; 35 36 ierr = DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);CHKERRQ(ierr); 37 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject) facegeom);CHKERRQ(ierr); 38 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject) cellgeom);CHKERRQ(ierr); 39 ierr = VecDestroy(&facegeom);CHKERRQ(ierr); 40 ierr = VecDestroy(&cellgeom);CHKERRQ(ierr); 41 } 42 if (facegeom) {PetscValidPointer(facegeom, 2); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject *) facegeom);CHKERRQ(ierr);} 43 if (cellgeom) {PetscValidPointer(cellgeom, 3); ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject *) cellgeom);CHKERRQ(ierr);} 44 if (minRadius) {ierr = DMPlexGetMinRadius(dm, minRadius);CHKERRQ(ierr);} 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "DMPlexTSGetGradientDM" 50 /*@C 51 DMPlexTSGetGradientDM - Return gradient data layout 52 53 Input Parameters: 54 + dm - The DM 55 - fv - The PetscFV 56 57 Output Parameter: 58 . dmGrad - The layout for gradient values 59 60 Level: developer 61 62 .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal() 63 @*/ 64 PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad) 65 { 66 DMTS dmts; 67 PetscObject obj; 68 PetscBool computeGradients; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 73 PetscValidHeaderSpecific(fv,PETSCFV_CLASSID,2); 74 PetscValidPointer(dmGrad,3); 75 ierr = PetscFVGetComputeGradients(fv, &computeGradients);CHKERRQ(ierr); 76 if (!computeGradients) {*dmGrad = NULL; PetscFunctionReturn(0);} 77 ierr = DMGetDMTS(dm, &dmts);CHKERRQ(ierr); 78 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);CHKERRQ(ierr); 79 if (!obj) { 80 DM dmGrad; 81 Vec faceGeometry, cellGeometry; 82 83 ierr = DMPlexTSGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);CHKERRQ(ierr); 84 ierr = DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);CHKERRQ(ierr); 85 ierr = PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);CHKERRQ(ierr); 86 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 87 } 88 ierr = PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "DMPlexTSComputeRHSFunctionFVM" 94 /*@ 95 DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user 96 97 Input Parameters: 98 + dm - The mesh 99 . t - The time 100 . locX - Local solution 101 - user - The user context 102 103 Output Parameter: 104 . F - Global output vector 105 106 Level: developer 107 108 .seealso: DMPlexComputeJacobianActionFEM() 109 @*/ 110 PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user) 111 { 112 Vec locF; 113 PetscInt cStart, cEnd, cEndInterior; 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 118 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 119 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 120 ierr = DMGetLocalVector(dm, &locF);CHKERRQ(ierr); 121 ierr = VecZeroEntries(locF);CHKERRQ(ierr); 122 ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, NULL, locF, user);CHKERRQ(ierr); 123 ierr = DMLocalToGlobalBegin(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr); 124 ierr = DMLocalToGlobalEnd(dm, locF, INSERT_VALUES, F);CHKERRQ(ierr); 125 ierr = DMRestoreLocalVector(dm, &locF);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "DMPlexTSComputeIFunctionFEM" 131 /*@ 132 DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user 133 134 Input Parameters: 135 + dm - The mesh 136 . t - The time 137 . locX - Local solution 138 . locX_t - Local solution time derivative, or NULL 139 - user - The user context 140 141 Output Parameter: 142 . locF - Local output vector 143 144 Level: developer 145 146 .seealso: DMPlexComputeJacobianActionFEM() 147 @*/ 148 PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user) 149 { 150 PetscInt cStart, cEnd, cEndInterior; 151 PetscErrorCode ierr; 152 153 PetscFunctionBegin; 154 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 155 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 156 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 157 ierr = DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, locX_t, locF, user);CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 #undef __FUNCT__ 162 #define __FUNCT__ "DMTSCheckFromOptions" 163 PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, void (**exactFuncs)(const PetscReal x[], PetscScalar *u, void *ctx), void **ctxs) 164 { 165 DM dm; 166 SNES snes; 167 Mat J, M; 168 Vec sol, r, b; 169 MatNullSpace nullSpace; 170 PetscReal *error, res = 0.0; 171 PetscInt numFields; 172 PetscBool check; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 ierr = PetscOptionsHasName(ts->hdr.prefix, "-dmts_check", &check);CHKERRQ(ierr); 177 if (!check) PetscFunctionReturn(0); 178 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 179 ierr = VecDuplicate(u, &sol);CHKERRQ(ierr); 180 ierr = VecDuplicate(u, &r);CHKERRQ(ierr); 181 ierr = TSSetSolution(ts, sol);CHKERRQ(ierr); 182 ierr = TSSetUp(ts);CHKERRQ(ierr); 183 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr); 184 ierr = SNESSetSolution(snes, sol);CHKERRQ(ierr); 185 ierr = DMCreateMatrix(dm, &J);CHKERRQ(ierr); 186 M = J; 187 #if 0 188 { 189 ierr = CreatePressureNullSpace(dm, user, &nullSpace);CHKERRQ(ierr); 190 ierr = MatSetNullSpace(J, nullSpace);CHKERRQ(ierr); 191 if (M != J) {ierr = MatSetNullSpace(M, nullSpace);CHKERRQ(ierr);} 192 } 193 #endif 194 /* Check discretization error */ 195 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 196 ierr = PetscMalloc1(PetscMax(1, numFields), &error);CHKERRQ(ierr); 197 if (numFields > 1) { 198 PetscInt f; 199 200 ierr = DMPlexComputeL2FieldDiff(dm, exactFuncs, ctxs, u, error);CHKERRQ(ierr); 201 ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");CHKERRQ(ierr); 202 for (f = 0; f < numFields; ++f) { 203 if (f) {ierr = PetscPrintf(PETSC_COMM_WORLD, ", ");CHKERRQ(ierr);} 204 if (error[f] >= 1.0e-11) {ierr = PetscPrintf(PETSC_COMM_WORLD, "%g", error[f]);CHKERRQ(ierr);} 205 else {ierr = PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");CHKERRQ(ierr);} 206 } 207 ierr = PetscPrintf(PETSC_COMM_WORLD, "]\n");CHKERRQ(ierr); 208 } else { 209 ierr = DMPlexComputeL2Diff(dm, exactFuncs, ctxs, u, &error[0]);CHKERRQ(ierr); 210 if (error[0] >= 1.0e-11) {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error[0]);CHKERRQ(ierr);} 211 else {ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");CHKERRQ(ierr);} 212 } 213 ierr = PetscFree(error);CHKERRQ(ierr); 214 /* Check residual */ 215 ierr = SNESComputeFunction(snes, u, r);CHKERRQ(ierr); 216 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 217 ierr = PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);CHKERRQ(ierr); 218 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr); 219 ierr = PetscObjectSetName((PetscObject) r, "Initial Residual");CHKERRQ(ierr); 220 ierr = VecViewFromOptions(r, "res_", "-vec_view");CHKERRQ(ierr); 221 /* Check Jacobian */ 222 ierr = SNESComputeJacobian(snes, u, M, M);CHKERRQ(ierr); 223 ierr = MatGetNullSpace(J, &nullSpace);CHKERRQ(ierr); 224 if (nullSpace) { 225 PetscBool isNull; 226 ierr = MatNullSpaceTest(nullSpace, J, &isNull);CHKERRQ(ierr); 227 if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid."); 228 } 229 ierr = VecDuplicate(u, &b);CHKERRQ(ierr); 230 ierr = VecSet(r, 0.0);CHKERRQ(ierr); 231 ierr = SNESComputeFunction(snes, r, b);CHKERRQ(ierr); 232 ierr = MatMult(M, u, r);CHKERRQ(ierr); 233 ierr = VecAXPY(r, 1.0, b);CHKERRQ(ierr); 234 ierr = VecDestroy(&b);CHKERRQ(ierr); 235 ierr = VecNorm(r, NORM_2, &res);CHKERRQ(ierr); 236 ierr = PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);CHKERRQ(ierr); 237 ierr = VecChop(r, 1.0e-10);CHKERRQ(ierr); 238 ierr = PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");CHKERRQ(ierr); 239 ierr = VecViewFromOptions(r, "linear_res_", "-vec_view");CHKERRQ(ierr); 240 ierr = VecDestroy(&sol);CHKERRQ(ierr); 241 ierr = VecDestroy(&r);CHKERRQ(ierr); 242 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 243 ierr = MatDestroy(&J);CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246