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