1 #include <petsc.h>
2
3 EXTERN_C_BEGIN
4 extern void formInitial(int*,int*,int*,double*,
5 double*,double*);
6 extern void formFunction(const int*,const int*,const int*,const double*,
7 const double*,const double[],const double[],double[]);
8 EXTERN_C_END
9
10 typedef struct AppCtx {
11 PetscInt nx,ny,nz;
12 PetscScalar h[3];
13 } AppCtx;
14
FormInitial(PetscReal t,Vec X,PetscCtx ctx)15 PetscErrorCode FormInitial(PetscReal t, Vec X, PetscCtx ctx)
16 {
17 PetscScalar *x;
18 AppCtx *app = (AppCtx*) ctx;
19
20 PetscFunctionBegin;
21 PetscCall(VecGetArray(X,&x));
22 /**/
23 formInitial(&app->nx,&app->ny,&app->nz,app->h,&t,x);
24 /**/
25 PetscCall(VecRestoreArray(X,&x));
26 PetscFunctionReturn(PETSC_SUCCESS);
27 }
28
FormFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)29 PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec Xdot,Vec F, PetscCtx ctx)
30 {
31 const PetscScalar *x;
32 const PetscScalar *xdot;
33 PetscScalar *f;
34 AppCtx *app = (AppCtx*) ctx;
35
36 PetscFunctionBegin;
37 PetscCall(VecGetArrayRead(X,&x));
38 PetscCall(VecGetArrayRead(Xdot,&xdot));
39 PetscCall(VecGetArray(F,&f));
40 /**/
41 formFunction(&app->nx,&app->ny,&app->nz,app->h,&t,x,xdot,f);
42 /**/
43 PetscCall(VecRestoreArrayRead(X,&x));
44 PetscCall(VecRestoreArrayRead(Xdot,&xdot));
45 PetscCall(VecRestoreArray(F,&f));
46 PetscFunctionReturn(PETSC_SUCCESS);
47 }
48
RunTest(int nx,int ny,int nz,int loops,double * wt)49 PetscErrorCode RunTest(int nx, int ny, int nz, int loops, double *wt)
50 {
51 Vec x,f;
52 TS ts;
53 AppCtx _app,*app=&_app;
54 double t1,t2;
55
56 PetscFunctionBegin;
57 app->nx = nx; app->h[0] = 1./(nx-1);
58 app->ny = ny; app->h[1] = 1./(ny-1);
59 app->nz = nz; app->h[2] = 1./(nz-1);
60
61 PetscCall(VecCreate(PETSC_COMM_SELF,&x));
62 PetscCall(VecSetSizes(x,nx*ny*nz,nx*ny*nz));
63 PetscCall(VecSetUp(x));
64 PetscCall(VecDuplicate(x,&f));
65
66 PetscCall(TSCreate(PETSC_COMM_SELF,&ts));
67 PetscCall(TSSetProblemType(ts,TS_NONLINEAR));
68 PetscCall(TSSetType(ts,TSTHETA));
69 PetscCall(TSThetaSetTheta(ts,1.0));
70 PetscCall(TSSetTimeStep(ts,0.01));
71 PetscCall(TSSetTime(ts,0.0));
72 PetscCall(TSSetMaxTime(ts,1.0));
73 PetscCall(TSSetMaxSteps(ts,10));
74 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER));
75
76 PetscCall(TSSetSolution(ts,x));
77 PetscCall(TSSetIFunction(ts,f,FormFunction,app));
78 PetscCall(PetscOptionsSetValue(NULL,"-snes_mf","1"));
79 {
80 SNES snes;
81 KSP ksp;
82 PetscCall(TSGetSNES(ts,&snes));
83 PetscCall(SNESGetKSP(snes,&ksp));
84 PetscCall(KSPSetType(ksp,KSPCG));
85 }
86 PetscCall(TSSetFromOptions(ts));
87 PetscCall(TSSetUp(ts));
88
89 *wt = 1e300;
90 while (loops-- > 0) {
91 PetscCall(FormInitial(0.0,x,app));
92 PetscCall(PetscTime(&t1));
93 PetscCall(TSSolve(ts,x));
94 PetscCall(PetscTime(&t2));
95 *wt = PetscMin(*wt,t2-t1);
96 }
97
98 PetscCall(VecDestroy(&x));
99 PetscCall(VecDestroy(&f));
100 PetscCall(TSDestroy(&ts));
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
GetInt(const char * name,PetscInt * v,PetscInt defv)104 PetscErrorCode GetInt(const char* name, PetscInt *v, PetscInt defv)
105 {
106 PetscFunctionBegin;
107 *v = defv;
108 PetscCall(PetscOptionsGetInt(NULL,NULL,name,v,NULL));
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111
main(int argc,char * argv[])112 int main(int argc, char *argv[])
113 {
114 double wt;
115 PetscInt n,start,step,stop,samples;
116
117 PetscCall(PetscInitialize(&argc,&argv,NULL,NULL));
118
119 PetscCall(GetInt("-start", &start, 12));
120 PetscCall(GetInt("-step", &step, 4));
121 PetscCall(GetInt("-stop", &stop, start));
122 PetscCall(GetInt("-samples", &samples, 1));
123
124 for (n=start; n<=stop; n+=step) {
125 int nx=n+1, ny=n+1, nz=n+1;
126 PetscCall(RunTest(nx,ny,nz,samples,&wt));
127 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Grid %3d x %3d x %3d -> %f seconds (%2d samples)\n",nx,ny,nz,wt,samples));
128 }
129 PetscCall(PetscFinalize());
130 return 0;
131 }
132