1 static char help[] = "Basic problem for multi-rate method.\n"; 2 3 /*F 4 5 \begin{eqnarray} 6 ys' = -2.0*\frac{-1.0+ys^2.0-\cos(t)}{2.0*ys}+0.05*\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-\frac{\sin(t)}{2.0*ys}\\ 7 yf' = 0.05*\frac{-1.0+ys^2-\cos(t)}{2.0*ys}-\frac{-2.0+yf^2-\cos(5.0*t)}{2.0*yf}-5.0*\frac{\sin(5.0*t)}{2.0*yf}\\ 8 \end{eqnarray} 9 10 F*/ 11 12 #include <petscts.h> 13 14 typedef struct { 15 PetscReal Tf,dt; 16 } AppCtx; 17 18 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 19 { 20 const PetscScalar *u; 21 PetscScalar *f; 22 23 PetscFunctionBegin; 24 PetscCall(VecGetArrayRead(U,&u)); 25 PetscCall(VecGetArray(F,&f)); 26 f[0] = -2.0*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])+0.05*(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-PetscSinScalar(t)/(2.0*u[0]); 27 f[1] = 0.05*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])-(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-5.0*PetscSinScalar(5.0*t)/(2.0*u[1]); 28 PetscCall(VecRestoreArrayRead(U,&u)); 29 PetscCall(VecRestoreArray(F,&f)); 30 PetscFunctionReturn(0); 31 } 32 33 static PetscErrorCode RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 34 { 35 const PetscScalar *u; 36 PetscScalar *f; 37 38 PetscFunctionBegin; 39 PetscCall(VecGetArrayRead(U,&u)); 40 PetscCall(VecGetArray(F,&f)); 41 f[0] = -2.0*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])+0.05*(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-PetscSinScalar(t)/(2.0*u[0]); 42 PetscCall(VecRestoreArrayRead(U,&u)); 43 PetscCall(VecRestoreArray(F,&f)); 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx) 48 { 49 const PetscScalar *u; 50 PetscScalar *f; 51 52 PetscFunctionBegin; 53 PetscCall(VecGetArrayRead(U,&u)); 54 PetscCall(VecGetArray(F,&f)); 55 f[0] = 0.05*(-1.0+u[0]*u[0]-PetscCosScalar(t))/(2.0*u[0])-(-2.0+u[1]*u[1]-PetscCosScalar(5.0*t))/(2.0*u[1])-5.0*PetscSinScalar(5.0*t)/(2.0*u[1]); 56 PetscCall(VecRestoreArrayRead(U,&u)); 57 PetscCall(VecRestoreArray(F,&f)); 58 PetscFunctionReturn(0); 59 } 60 61 static PetscErrorCode sol_true(PetscReal t,Vec U) 62 { 63 PetscScalar *u; 64 65 PetscFunctionBegin; 66 PetscCall(VecGetArray(U,&u)); 67 u[0] = PetscSqrtScalar(1.0+PetscCosScalar(t)); 68 u[1] = PetscSqrtScalar(2.0+PetscCosScalar(5.0*t)); 69 PetscCall(VecRestoreArray(U,&u)); 70 PetscFunctionReturn(0); 71 } 72 73 int main(int argc,char **argv) 74 { 75 TS ts; /* ODE integrator */ 76 Vec U; /* solution will be stored here */ 77 Vec Utrue; 78 PetscMPIInt size; 79 AppCtx ctx; 80 PetscScalar *u; 81 IS iss; 82 IS isf; 83 PetscInt *indicess; 84 PetscInt *indicesf; 85 PetscInt n=2; 86 PetscReal error,tt; 87 88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 89 Initialize program 90 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 91 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 92 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 93 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Create index for slow part and fast part 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 PetscCall(PetscMalloc1(1,&indicess)); 99 indicess[0]=0; 100 PetscCall(PetscMalloc1(1,&indicesf)); 101 indicesf[0]=1; 102 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss)); 103 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf)); 104 105 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 106 Create necessary vector 107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 108 PetscCall(VecCreate(PETSC_COMM_WORLD,&U)); 109 PetscCall(VecSetSizes(U,n,PETSC_DETERMINE)); 110 PetscCall(VecSetFromOptions(U)); 111 PetscCall(VecDuplicate(U,&Utrue)); 112 PetscCall(VecCopy(U,Utrue)); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Set initial condition 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 PetscCall(VecGetArray(U,&u)); 118 u[0] = PetscSqrtScalar(2.0); 119 u[1] = PetscSqrtScalar(3.0); 120 PetscCall(VecRestoreArray(U,&u)); 121 122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 123 Create timestepping solver context 124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 125 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 126 PetscCall(TSSetType(ts,TSMPRK)); 127 128 PetscCall(TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx)); 129 PetscCall(TSRHSSplitSetIS(ts,"slow",iss)); 130 PetscCall(TSRHSSplitSetIS(ts,"fast",isf)); 131 PetscCall(TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx)); 132 PetscCall(TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx)); 133 134 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 135 Set initial conditions 136 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 137 PetscCall(TSSetSolution(ts,U)); 138 139 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 140 Set solver options 141 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 142 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options",""); 143 { 144 ctx.Tf = 0.3; 145 ctx.dt = 0.01; 146 PetscCall(PetscOptionsScalar("-Tf","","",ctx.Tf,&ctx.Tf,NULL)); 147 PetscCall(PetscOptionsScalar("-dt","","",ctx.dt,&ctx.dt,NULL)); 148 } 149 PetscOptionsEnd(); 150 PetscCall(TSSetMaxTime(ts,ctx.Tf)); 151 PetscCall(TSSetTimeStep(ts,ctx.dt)); 152 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 153 PetscCall(TSSetFromOptions(ts)); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 Solve linear system 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 PetscCall(TSSolve(ts,U)); 159 PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); 160 161 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162 Check the error of the Petsc solution 163 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 164 PetscCall(TSGetTime(ts,&tt)); 165 PetscCall(sol_true(tt,Utrue)); 166 PetscCall(VecAXPY(Utrue,-1.0,U)); 167 PetscCall(VecNorm(Utrue,NORM_2,&error)); 168 169 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 170 Print norm2 error 171 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 172 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"l2 error norm: %g\n", (double)error)); 173 174 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 175 Free work space. All PETSc objects should be destroyed when they are no longer needed. 176 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 177 PetscCall(VecDestroy(&U)); 178 PetscCall(TSDestroy(&ts)); 179 PetscCall(VecDestroy(&Utrue)); 180 PetscCall(ISDestroy(&iss)); 181 PetscCall(ISDestroy(&isf)); 182 PetscCall(PetscFree(indicess)); 183 PetscCall(PetscFree(indicesf)); 184 PetscCall(PetscFinalize()); 185 return 0; 186 } 187 188 /*TEST 189 build: 190 requires: !complex 191 192 test: 193 194 TEST*/ 195