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
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)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
RHSFunctionslow(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)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
RHSFunctionfast(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)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
sol_true(PetscReal t,Vec U)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
main(int argc,char ** argv)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
196
197 test:
198
199 TEST*/
200