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,(char*)0,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: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[ctxt->physics_type]); 137 } 138 } 139 PetscCall(VecRestoreArray(u,&uarr)); 140 PetscFunctionReturn(0); 141 } 142 143 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat Jpre,void *ctx) 144 { 145 UserContext *user = (UserContext*) ctx; 146 PetscInt matis,matie,i; 147 PetscReal dx,dx2; 148 149 PetscFunctionBeginUser; 150 dx = (user->xmax - user->xmin)/((PetscReal)user->imax); dx2 = dx*dx; 151 PetscCall(MatGetOwnershipRange(J,&matis,&matie)); 152 for (i=matis; i<matie; i++) { 153 PetscScalar values[3]; 154 PetscInt col[3]; 155 switch (user->physics_type) { 156 case PHYSICS_DIFFUSION: 157 values[0] = user->a*1.0/dx2; 158 values[1] = -user->a*2.0/dx2; 159 values[2] = user->a*1.0/dx2; 160 break; 161 case PHYSICS_ADVECTION: 162 values[0] = user->a*.5/dx; 163 values[1] = 0.; 164 values[2] = -user->a*.5/dx; 165 break; 166 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[user->physics_type]); 167 } 168 /* periodic boundaries */ 169 if (i == 0) { 170 col[0] = user->imax-1; 171 col[1] = i; 172 col[2] = i+1; 173 } else if (i == user->imax-1) { 174 col[0] = i-1; 175 col[1] = i; 176 col[2] = 0; 177 } else { 178 col[0] = i-1; 179 col[1] = i; 180 col[2] = i+1; 181 } 182 PetscCall(MatSetValues(J,1,&i,3,col,values,INSERT_VALUES)); 183 } 184 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 185 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 186 PetscFunctionReturn(0); 187 } 188 189 /*TEST 190 191 test: 192 requires: double 193 suffix: 1 194 nsize: {{1 2}} 195 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 2 196 197 test: 198 requires: double 199 suffix: 2 200 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 3 201 202 testset: 203 requires: hpddm 204 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 205 test: 206 suffix: 3 207 requires: double 208 args: -ksp_hpddm_precision {{single double}shared output} 209 test: 210 suffix: 3_single 211 requires: single 212 args: -ksp_hpddm_precision {{single double}shared output} 213 test: 214 suffix: 3___float128 215 requires: __float128 216 output_file: output/ex74_3.out 217 args: -ksp_hpddm_precision {{double quadruple}shared output} 218 TEST*/ 219