Actual source code: ex1.c
petsc-3.13.4 2020-08-01
1: static char help[] = "Basic problem for multi-rate method.\n";
\begin{eqnarray}
ys' = \frac{ys}{a}\\
yf' = ys*cos(bt)\\
\end{eqnarray}
12: #include <petscts.h>
14: typedef struct {
15: PetscReal a,b,Tf,dt;
16: } AppCtx;
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;
25: VecGetArrayRead(U,&u);
26: VecGetArray(F,&f);
27: f[0] = u[0]/ctx->a;
28: f[1] = u[0]*PetscCosScalar(t*ctx->b);
29: VecRestoreArrayRead(U,&u);
30: VecRestoreArray(F,&f);
31: return(0);
32: }
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;
41: VecGetArrayRead(U,&u);
42: VecGetArray(F,&f);
43: f[0] = u[0]/ctx->a;
44: VecRestoreArrayRead(U,&u);
45: VecRestoreArray(F,&f);
46: return(0);
47: }
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;
56: VecGetArrayRead(U,&u);
57: VecGetArray(F,&f);
58: f[0] = u[0]*PetscCosScalar(t*ctx->b);
59: VecRestoreArrayRead(U,&u);
60: VecRestoreArray(F,&f);
61: return(0);
62: }
64: /*
65: Define the analytic solution for check method easily
66: */
67: static PetscErrorCode sol_true(PetscReal t, Vec U,AppCtx *ctx)
68: {
70: PetscScalar *u;
73: VecGetArray(U,&u);
74: u[0] = PetscExpScalar(t/ctx->a);
75: u[1] = (ctx->a*PetscCosScalar(ctx->b*t)+ctx->a*ctx->a*ctx->b*PetscSinScalar(ctx->b*t))*PetscExpScalar(t/ctx->a)/(1+ctx->a*ctx->a*ctx->b*ctx->b);
76: VecRestoreArray(U,&u);
77: return(0);
78: }
80: int main(int argc,char **argv)
81: {
82: TS ts; /* ODE integrator */
83: Vec U; /* solution will be stored here */
84: Vec Utrue;
86: PetscMPIInt size;
87: AppCtx ctx;
88: PetscScalar *u;
89: IS iss;
90: IS isf;
91: PetscInt *indicess;
92: PetscInt *indicesf;
93: PetscInt n=2;
94: PetscReal error,tt;
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Initialize program
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
100: MPI_Comm_size(PETSC_COMM_WORLD,&size);
101: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create index for slow part and fast part
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: PetscMalloc1(1,&indicess);
107: indicess[0] = 0;
108: PetscMalloc1(1,&indicesf);
109: indicesf[0] = 1;
110: ISCreateGeneral(PETSC_COMM_SELF,1,indicess,PETSC_COPY_VALUES,&iss);
111: ISCreateGeneral(PETSC_COMM_SELF,1,indicesf,PETSC_COPY_VALUES,&isf);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Create necesary vector
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: VecCreate(PETSC_COMM_WORLD,&U);
117: VecSetSizes(U,n,PETSC_DETERMINE);
118: VecSetFromOptions(U);
119: VecDuplicate(U,&Utrue);
120: VecCopy(U,Utrue);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Set runtime options
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ODE options","");
126: {
127: ctx.a = 2.0;
128: ctx.b = 25.0;
129: PetscOptionsReal("-a","","",ctx.a,&ctx.a,NULL);
130: PetscOptionsReal("-b","","",ctx.b,&ctx.b,NULL);
131: ctx.Tf = 2;
132: ctx.dt = 0.01;
133: PetscOptionsReal("-Tf","","",ctx.Tf,&ctx.Tf,NULL);
134: PetscOptionsReal("-dt","","",ctx.dt,&ctx.dt,NULL);
135: }
136: PetscOptionsEnd();
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Initialize the solution
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: VecGetArray(U,&u);
142: u[0] = 1.0;
143: u[1] = ctx.a/(1+ctx.a*ctx.a*ctx.b*ctx.b);
144: VecRestoreArray(U,&u);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Create timestepping solver context
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: TSCreate(PETSC_COMM_WORLD,&ts);
150: TSSetType(ts,TSMPRK);
152: TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
153: TSRHSSplitSetIS(ts,"slow",iss);
154: TSRHSSplitSetIS(ts,"fast",isf);
155: TSRHSSplitSetRHSFunction(ts,"slow",NULL,(TSRHSFunction)RHSFunctionslow,&ctx);
156: TSRHSSplitSetRHSFunction(ts,"fast",NULL,(TSRHSFunction)RHSFunctionfast,&ctx);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Set initial conditions
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: TSSetSolution(ts,U);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Set solver options
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: TSSetMaxTime(ts,ctx.Tf);
167: TSSetTimeStep(ts,ctx.dt);
168: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
169: TSSetFromOptions(ts);
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Solve linear system
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: TSSolve(ts,U);
175: VecView(U,PETSC_VIEWER_STDOUT_WORLD);
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Check the error of the Petsc solution
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: TSGetTime(ts,&tt);
181: sol_true(tt,Utrue,&ctx);
182: VecAXPY(Utrue,-1.0,U);
183: VecNorm(Utrue,NORM_2,&error);
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Print norm2 error
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188: PetscPrintf(PETSC_COMM_WORLD,"l2 error norm = %g\n", error);
190: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191: Free work space. All PETSc objects should be destroyed when they are no longer needed.
192: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193: VecDestroy(&U);
194: VecDestroy(&Utrue);
195: TSDestroy(&ts);
196: ISDestroy(&iss);
197: ISDestroy(&isf);
198: PetscFree(indicess);
199: PetscFree(indicesf);
200: PetscFinalize();
201: return ierr;
202: }
204: /*TEST
205: build:
206: requires: !complex c99
208: test:
210: TEST*/