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 PetscErrorCode ierr; 53 54 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 55 56 /* default value */ 57 ctxt.a = 0.1; 58 ctxt.xmin = 0.0; 59 ctxt.xmax = 1.0; 60 ctxt.imax = 40; 61 ctxt.physics_type = PHYSICS_DIFFUSION; 62 63 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"IRK options","");CHKERRQ(ierr); 64 ierr = PetscOptionsReal("-a","diffusion coefficient","<1.0>",ctxt.a,&ctxt.a,NULL);CHKERRQ(ierr); 65 ierr = PetscOptionsInt ("-imax","grid size","<20>",ctxt.imax,&ctxt.imax,NULL);CHKERRQ(ierr); 66 ierr = PetscOptionsReal("-xmin","xmin","<0.0>",ctxt.xmin,&ctxt.xmin,NULL);CHKERRQ(ierr); 67 ierr = PetscOptionsReal("-xmax","xmax","<1.0>",ctxt.xmax,&ctxt.xmax,NULL);CHKERRQ(ierr); 68 ierr = PetscOptionsEnum("-physics_type","Type of process to discretize","",PhysicsTypes,(PetscEnum)ctxt.physics_type,(PetscEnum*)&ctxt.physics_type,NULL);CHKERRQ(ierr); 69 ierr = PetscOptionsEnd();CHKERRQ(ierr); 70 71 /* allocate and initialize solution vector and exact solution */ 72 ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr); 73 ierr = VecSetSizes(u,PETSC_DECIDE,ctxt.imax);CHKERRQ(ierr); 74 ierr = VecSetFromOptions(u);CHKERRQ(ierr); 75 ierr = VecDuplicate(u,&uex);CHKERRQ(ierr); 76 /* initial solution */ 77 ierr = ExactSolution(u,&ctxt,0.0);CHKERRQ(ierr); 78 79 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 80 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 81 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,ctxt.imax,ctxt.imax);CHKERRQ(ierr); 82 ierr = MatSetUp(A);CHKERRQ(ierr); 83 84 /* Create and set options for TS */ 85 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 86 ierr = TSSetProblemType(ts,TS_LINEAR);CHKERRQ(ierr); 87 ierr = TSSetTimeStep(ts,0.125);CHKERRQ(ierr); 88 ierr = TSSetSolution(ts,u);CHKERRQ(ierr); 89 ierr = TSSetMaxSteps(ts,10);CHKERRQ(ierr); 90 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr); 91 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 92 ierr = TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&ctxt);CHKERRQ(ierr); 93 ierr = RHSJacobian(ts,0,u,A,A,&ctxt);CHKERRQ(ierr); 94 ierr = TSSetRHSJacobian(ts,A,A,TSComputeRHSJacobianConstant,&ctxt);CHKERRQ(ierr); 95 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 96 ierr = TSSolve(ts,u);CHKERRQ(ierr); 97 98 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 99 /* exact solution */ 100 ierr = ExactSolution(uex,&ctxt,ftime);CHKERRQ(ierr); 101 102 /* Calculate error in final solution */ 103 ierr = VecAYPX(uex,-1.0,u);CHKERRQ(ierr); 104 ierr = VecNorm(uex,NORM_2,&err);CHKERRQ(ierr); 105 err = PetscSqrtReal(err*err/((PetscReal)ctxt.imax)); 106 ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 norm of the numerical error = %g (time=%g)\n",(double)err,(double)ftime);CHKERRQ(ierr); 107 108 /* Free up memory */ 109 ierr = TSDestroy(&ts);CHKERRQ(ierr); 110 ierr = MatDestroy(&A);CHKERRQ(ierr); 111 ierr = VecDestroy(&uex);CHKERRQ(ierr); 112 ierr = VecDestroy(&u);CHKERRQ(ierr); 113 ierr = PetscFinalize(); 114 return ierr; 115 } 116 117 PetscErrorCode ExactSolution(Vec u,void *c,PetscReal t) 118 { 119 UserContext *ctxt = (UserContext*) c; 120 PetscErrorCode ierr; 121 PetscInt i,is,ie; 122 PetscScalar *uarr; 123 PetscReal x,dx,a=ctxt->a,pi=PETSC_PI; 124 125 PetscFunctionBegin; 126 dx = (ctxt->xmax - ctxt->xmin)/((PetscReal) ctxt->imax); 127 ierr = VecGetOwnershipRange(u,&is,&ie);CHKERRQ(ierr); 128 ierr = VecGetArray(u,&uarr);CHKERRQ(ierr); 129 for (i=is; i<ie; i++) { 130 x = i * dx; 131 switch (ctxt->physics_type) { 132 case PHYSICS_DIFFUSION: 133 uarr[i-is] = PetscExpScalar(-4.0*pi*pi*a*t)*PetscSinScalar(2*pi*x); 134 break; 135 case PHYSICS_ADVECTION: 136 uarr[i-is] = PetscSinScalar(2*pi*(x - a*t)); 137 break; 138 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[ctxt->physics_type]); 139 } 140 } 141 ierr = VecRestoreArray(u,&uarr);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat Jpre,void *ctx) 146 { 147 UserContext *user = (UserContext*) ctx; 148 PetscInt matis,matie,i; 149 PetscReal dx,dx2; 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 dx = (user->xmax - user->xmin)/((PetscReal)user->imax); dx2 = dx*dx; 154 ierr = MatGetOwnershipRange(J,&matis,&matie);CHKERRQ(ierr); 155 for (i=matis; i<matie; i++) { 156 PetscScalar values[3]; 157 PetscInt col[3]; 158 switch (user->physics_type) { 159 case PHYSICS_DIFFUSION: 160 values[0] = user->a*1.0/dx2; 161 values[1] = -user->a*2.0/dx2; 162 values[2] = user->a*1.0/dx2; 163 break; 164 case PHYSICS_ADVECTION: 165 values[0] = user->a*.5/dx; 166 values[1] = 0.; 167 values[2] = -user->a*.5/dx; 168 break; 169 default: SETERRQ1(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 ierr = MatSetValues(J,1,&i,3,col,values,INSERT_VALUES);CHKERRQ(ierr); 186 } 187 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 188 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 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 test: 206 requires: hpddm double 207 suffix: 3 208 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 3 -ksp_view_final_residual -ksp_hpddm_type bgcrodr -ksp_type hpddm 209 TEST*/ 210