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