1 static char help[] = "Solves the constant-coefficient 1D heat equation \n\ 2 with an Implicit Runge-Kutta method using MatKAIJ. \n\ 3 \n\ 4 du d^2 u \n\ 5 -- = a ----- ; 0 <= x <= 1; \n\ 6 dt dx^2 \n\ 7 \n\ 8 with periodic boundary conditions \n\ 9 \n\ 10 2nd order central discretization in space: \n\ 11 \n\ 12 [ d^2 u ] u_{i+1} - 2u_i + u_{i-1} \n\ 13 [ ----- ] = ------------------------ \n\ 14 [ dx^2 ]i h^2 \n\ 15 \n\ 16 i = grid index; h = x_{i+1}-x_i (Uniform) \n\ 17 0 <= i < n h = 1.0/n \n\ 18 \n\ 19 Thus, \n\ 20 \n\ 21 du \n\ 22 -- = Ju; J = (a/h^2) tridiagonal(1,-2,1)_n \n\ 23 dt \n\ 24 \n\ 25 This example is a TS version of the KSP ex74.c tutorial. \n"; 26 27 #include <petscts.h> 28 29 typedef enum { 30 PHYSICS_DIFFUSION, 31 PHYSICS_ADVECTION 32 } PhysicsType; 33 const char *const PhysicsTypes[] = {"DIFFUSION", "ADVECTION", "PhysicsType", "PHYSICS_", NULL}; 34 35 typedef struct Context { 36 PetscReal a; /* diffusion coefficient */ 37 PetscReal xmin, xmax; /* domain bounds */ 38 PetscInt imax; /* number of grid points */ 39 PhysicsType physics_type; 40 } UserContext; 41 42 static PetscErrorCode ExactSolution(Vec, void *, PetscReal); 43 static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *); 44 45 int main(int argc, char **argv) 46 { 47 TS ts; 48 Mat A; 49 Vec u, uex; 50 UserContext ctxt; 51 PetscReal err, ftime; 52 53 PetscFunctionBeginUser; 54 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 55 /* default value */ 56 ctxt.a = 0.1; 57 ctxt.xmin = 0.0; 58 ctxt.xmax = 1.0; 59 ctxt.imax = 40; 60 ctxt.physics_type = PHYSICS_DIFFUSION; 61 62 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "IRK options", ""); 63 PetscCall(PetscOptionsReal("-a", "diffusion coefficient", "<1.0>", ctxt.a, &ctxt.a, NULL)); 64 PetscCall(PetscOptionsInt("-imax", "grid size", "<20>", ctxt.imax, &ctxt.imax, NULL)); 65 PetscCall(PetscOptionsReal("-xmin", "xmin", "<0.0>", ctxt.xmin, &ctxt.xmin, NULL)); 66 PetscCall(PetscOptionsReal("-xmax", "xmax", "<1.0>", ctxt.xmax, &ctxt.xmax, NULL)); 67 PetscCall(PetscOptionsEnum("-physics_type", "Type of process to discretize", "", PhysicsTypes, (PetscEnum)ctxt.physics_type, (PetscEnum *)&ctxt.physics_type, NULL)); 68 PetscOptionsEnd(); 69 70 /* allocate and initialize solution vector and exact solution */ 71 PetscCall(VecCreate(PETSC_COMM_WORLD, &u)); 72 PetscCall(VecSetSizes(u, PETSC_DECIDE, ctxt.imax)); 73 PetscCall(VecSetFromOptions(u)); 74 PetscCall(VecDuplicate(u, &uex)); 75 /* initial solution */ 76 PetscCall(ExactSolution(u, &ctxt, 0.0)); 77 78 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 79 PetscCall(MatSetType(A, MATAIJ)); 80 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, ctxt.imax, ctxt.imax)); 81 PetscCall(MatSetUp(A)); 82 83 /* Create and set options for TS */ 84 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 85 PetscCall(TSSetProblemType(ts, TS_LINEAR)); 86 PetscCall(TSSetTimeStep(ts, 0.125)); 87 PetscCall(TSSetSolution(ts, u)); 88 PetscCall(TSSetMaxSteps(ts, 10)); 89 PetscCall(TSSetMaxTime(ts, 1.0)); 90 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 91 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &ctxt)); 92 PetscCall(RHSJacobian(ts, 0, u, A, A, &ctxt)); 93 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &ctxt)); 94 PetscCall(TSSetFromOptions(ts)); 95 PetscCall(TSSolve(ts, u)); 96 97 PetscCall(TSGetSolveTime(ts, &ftime)); 98 /* exact solution */ 99 PetscCall(ExactSolution(uex, &ctxt, ftime)); 100 101 /* Calculate error in final solution */ 102 PetscCall(VecAYPX(uex, -1.0, u)); 103 PetscCall(VecNorm(uex, NORM_2, &err)); 104 err = PetscSqrtReal(err * err / ((PetscReal)ctxt.imax)); 105 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L2 norm of the numerical error = %g (time=%g)\n", (double)err, (double)ftime)); 106 107 /* Free up memory */ 108 PetscCall(TSDestroy(&ts)); 109 PetscCall(MatDestroy(&A)); 110 PetscCall(VecDestroy(&uex)); 111 PetscCall(VecDestroy(&u)); 112 PetscCall(PetscFinalize()); 113 return 0; 114 } 115 116 PetscErrorCode ExactSolution(Vec u, void *c, PetscReal t) 117 { 118 UserContext *ctxt = (UserContext *)c; 119 PetscInt i, is, ie; 120 PetscScalar *uarr; 121 PetscReal x, dx, a = ctxt->a, pi = PETSC_PI; 122 123 PetscFunctionBeginUser; 124 dx = (ctxt->xmax - ctxt->xmin) / ((PetscReal)ctxt->imax); 125 PetscCall(VecGetOwnershipRange(u, &is, &ie)); 126 PetscCall(VecGetArray(u, &uarr)); 127 for (i = is; i < ie; i++) { 128 x = i * dx; 129 switch (ctxt->physics_type) { 130 case PHYSICS_DIFFUSION: 131 uarr[i - is] = PetscExpScalar(-4.0 * pi * pi * a * t) * PetscSinScalar(2 * pi * x); 132 break; 133 case PHYSICS_ADVECTION: 134 uarr[i - is] = PetscSinScalar(2 * pi * (x - a * t)); 135 break; 136 default: 137 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for physics type %s", PhysicsTypes[ctxt->physics_type]); 138 } 139 } 140 PetscCall(VecRestoreArray(u, &uarr)); 141 PetscFunctionReturn(PETSC_SUCCESS); 142 } 143 144 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx) 145 { 146 UserContext *user = (UserContext *)ctx; 147 PetscInt matis, matie, i; 148 PetscReal dx, dx2; 149 150 PetscFunctionBeginUser; 151 dx = (user->xmax - user->xmin) / ((PetscReal)user->imax); 152 dx2 = dx * dx; 153 PetscCall(MatGetOwnershipRange(J, &matis, &matie)); 154 for (i = matis; i < matie; i++) { 155 PetscScalar values[3]; 156 PetscInt col[3]; 157 switch (user->physics_type) { 158 case PHYSICS_DIFFUSION: 159 values[0] = user->a * 1.0 / dx2; 160 values[1] = -user->a * 2.0 / dx2; 161 values[2] = user->a * 1.0 / dx2; 162 break; 163 case PHYSICS_ADVECTION: 164 values[0] = user->a * .5 / dx; 165 values[1] = 0.; 166 values[2] = -user->a * .5 / dx; 167 break; 168 default: 169 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for physics type %s", PhysicsTypes[user->physics_type]); 170 } 171 /* periodic boundaries */ 172 if (i == 0) { 173 col[0] = user->imax - 1; 174 col[1] = i; 175 col[2] = i + 1; 176 } else if (i == user->imax - 1) { 177 col[0] = i - 1; 178 col[1] = i; 179 col[2] = 0; 180 } else { 181 col[0] = i - 1; 182 col[1] = i; 183 col[2] = i + 1; 184 } 185 PetscCall(MatSetValues(J, 1, &i, 3, col, values, INSERT_VALUES)); 186 } 187 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 188 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 189 PetscFunctionReturn(PETSC_SUCCESS); 190 } 191 192 /*TEST 193 194 test: 195 requires: double 196 suffix: 1 197 nsize: {{1 2}} 198 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 2 199 200 test: 201 requires: double 202 suffix: 2 203 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 3 204 205 testset: 206 requires: hpddm 207 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-4 -ts_type irk -ts_irk_nstages 3 -ksp_view_final_residual -ksp_hpddm_type gcrodr -ksp_type hpddm 208 test: 209 suffix: 3 210 requires: double 211 args: -ksp_hpddm_precision {{single double}shared output} 212 test: 213 suffix: 3_single 214 requires: single 215 args: -ksp_hpddm_precision {{single double}shared output} 216 test: 217 suffix: 3___float128 218 requires: __float128 219 output_file: output/ex74_3.out 220 args: -ksp_hpddm_precision {{double quadruple}shared output} 221 TEST*/ 222