static char help[] = "Time-dependent 1D diffusion equation with FV discretization and cell-centered grid\n"; /* u_t = uxx, 0 < x < 1 with IC: u(x,t=0) = 1.0, if 0.4 < x < 0.6 u(x,t=0) = 0.0, otherwise and Dirichlet BC: u(x=0) = u_A = 0.0; u(x=1) = u_B = 1.0; or Neumann BC: du/dx(x=0.0) = q_A = 0.0; du/dx(x=1.0) = q_B = 1.0; mpirun -np 1 ./UnsteadyCond2D -da_grid_x 100 -ts_monitor_draw_solution -draw_pause 0.01 -ts_monitor */ #include #include #include /* User-defined data structures and routines*/ typedef struct { PetscInt nx, ny, nz; PetscReal L, M, N; PetscReal dx, dy, dz; PetscReal xmin, xmax, ymin, ymax, zmin, zmax; } AppCtx; extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void*); extern PetscErrorCode FormInitialSolution(DM, Vec, void*); int main(int argc,char **argv) { TS ts; /* nonlinear solver */ Vec u; /* solution vector */ PetscInt steps; /* iterations for convergence */ PetscErrorCode ierr; DM da; PetscReal ftime,dt; AppCtx user; /* user-defined work context */ DMDALocalInfo info; ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; user.xmin = 0.0; user.xmax = 1.0; user.ymin = 0.0; user.ymax = 1.0; user.zmin = 0.0; user.zmax = 0.0; user.nx = 20; // x-dir number of FV cells user.ny = 1; // y-dir number of FV cells user.nz = 1; // z-dir number of FV cells user.L = user.xmax - user.xmin; // x-dir size of domain user.M = user.ymax - user.ymin; // y-dir size of domain user.N = user.zmax - user.zmin; // z-dir size of domain ftime = 100.0; // max time dt = .01; // timestep ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.nx, 1, 1, NULL, &da); CHKERRQ(ierr); ierr = DMSetFromOptions(da); CHKERRQ(ierr); ierr = DMSetUp(da); CHKERRQ(ierr); ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr); // sets a DMDA coordinates to be a uniform cell-centered grid user.dx = user.L/info.mx; user.dy = user.M/info.my; user.dz = user.N/info.mz; ierr = PetscPrintf(PETSC_COMM_WORLD, "user.dx, dy, dz = %.3g, %.3g, %.3g\n", user.dx, user.dy, user.dz);CHKERRQ(ierr); // transfer a vertex-centered grid to cell-centered FV gird ierr = DMDASetUniformCoordinates(da, user.xmin + user.dx / 2.0, user.xmax - user.dx / 2.0, // x-dir extremes user.ymin + user.dy / 2.0, user.ymax - user.dy / 2.0, // y-dir extremes user.zmin + user.dz / 2.0, user.zmax - user.dz / 2.0 // z-dir extremes ); ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr); ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); ierr = TSSetDM(ts,da);CHKERRQ(ierr); ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr); //ierr = TSSetType(ts,TSBDF);CHKERRQ(ierr); ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr); ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); // initial condition ierr = FormInitialSolution(da,u,&user);CHKERRQ(ierr); ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); ierr = TSSetFromOptions(ts);CHKERRQ(ierr); ierr = TSSolve(ts,u);CHKERRQ(ierr); ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = TSDestroy(&ts);CHKERRQ(ierr); ierr = DMDestroy(&da);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr; } /* ------------------------------------------------------------------- */ /* RHSFunction - Evaluates nonlinear function, G(u). Input Parameters: . ts - the TS context . U - input vector . ptr - optional user-defined context, as set by TSSetFunction() Output Parameter: . G - function vector*/ PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec U, Vec G, void *ptr) { //AppCtx *user=(AppCtx*)ptr; DM da; PetscErrorCode ierr; PetscInt i, xs, xm, mx; PetscReal dx; PetscScalar *u, *f; Vec localU; DMDALocalInfo info; PetscReal u_A, u_B, u_C, u_D; //Dirichlet BC PetscReal q_A, q_B, q_C, q_D; //Neumann BC u_A = 0.0; u_B = 0.0; //Dirichlet BC q_A = 0.0; q_B = 0.0; //Neumann BC, q = du/dx PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr); ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr); xs = info.xs; xm = info.xm; mx = info.mx; dx = 1.0 / (PetscReal) (mx); ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU); CHKERRQ(ierr); ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU); CHKERRQ(ierr); ierr = DMDAVecGetArrayRead(da,localU,&u); CHKERRQ(ierr); ierr = DMDAVecGetArray(da,G,&f); CHKERRQ(ierr); for (i = xs; i < xs + xm; i++) { if (i == 0) { //Dirichlet BC f[i] = u[i]; //f[i] = u[i] - u_A; // dose not work //Neumann BC //f[i] = 3.0 * u[i] - u[i+1] - 2.0 * u_A; // dose not work //f[i] = (u[i]-u[i-1])/dx + q_A; // dose not work continue; } else if (i == mx - 1) { //Dirichlet BC f[i] = u[i]; //f[i] = u[i] - u_B; // not work //Neumann BC //f[i] = 3.0 * u[i] - u[i-1] - 2.0 * u_B; // dose not work //f[i] = (u[i-1]-u[i])/dx + q_B; // dose not work continue; } else { f[i] = (u[i - 1] - 2.0 * u[i] + u[i + 1]) / dx / dx; } } ierr = DMDAVecRestoreArrayRead(da,localU,&u); CHKERRQ(ierr); ierr = DMDAVecRestoreArray(da,G,&f); CHKERRQ(ierr); ierr = DMRestoreLocalVector(da,&localU); CHKERRQ(ierr); PetscFunctionReturn(0); } /* ------------------------------------------------------------------- */ /* A squre wave is set as the initial condition */ PetscErrorCode FormInitialSolution(DM da, Vec U, void *ptr) { AppCtx *user=(AppCtx*)ptr; PetscErrorCode ierr; PetscInt i,xs,xm,mx; PetscScalar *u; PetscReal dx,x; DMDALocalInfo info; PetscFunctionBeginUser; ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr); xs = info.xs; xm = info.xm; mx = info.mx; dx = 1.0 / (PetscReal)(mx); ierr = DMDAGetLocalInfo(da, &info); CHKERRQ(ierr); xs = info.xs; xm = info.xm; mx = info.mx; // u = 1.0, 0.4 < x < 0.6 // u = 0.0, otherwise ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); for (i = xs; i < xs + xm; i++) { x = dx / 2 + dx * i; if (x > 0.4 && x < 0.6) u[i] = 1.0; else u[i] = 0.0; } ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); PetscFunctionReturn(0); }