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 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 29 PetscErrorCode FormFunction(TS ts, PetscReal t, Vec X, Vec Xdot,Vec F, void *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 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 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 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