1 static char help[] = "Solves the time dependent Allen-Cahn equation with IMEX methods"; 2 3 /* 4 * allen_cahn.c 5 * 6 * Created on: Jun 8, 2012 7 * Author: Hong Zhang 8 */ 9 10 #include <petscts.h> 11 12 /* 13 * application context 14 */ 15 typedef struct { 16 PetscReal param; /* parameter */ 17 PetscReal xleft, xright; /* range in x-direction */ 18 PetscInt mx; /* Discretization in x-direction */ 19 } AppCtx; 20 21 static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 22 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *); 23 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *ctx); 24 static PetscErrorCode FormInitialSolution(TS, Vec, void *); 25 26 int main(int argc, char **argv) { 27 TS ts; 28 Vec x; /* solution vector */ 29 Mat A; /* Jacobian */ 30 PetscInt steps, mx; 31 PetscReal ftime; 32 AppCtx user; /* user-defined work context */ 33 34 PetscFunctionBeginUser; 35 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 36 37 /* Initialize user application context */ 38 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Allen-Cahn equation", ""); 39 user.param = 9e-4; 40 user.xleft = -1.; 41 user.xright = 2.; 42 user.mx = 400; 43 PetscCall(PetscOptionsReal("-eps", "parameter", "", user.param, &user.param, NULL)); 44 PetscOptionsEnd(); 45 46 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47 Set runtime options 48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 49 /* 50 * PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 51 */ 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Create necessary matrix and vectors, solve same ODE on every process 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 56 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, user.mx, user.mx)); 57 PetscCall(MatSetFromOptions(A)); 58 PetscCall(MatSetUp(A)); 59 PetscCall(MatCreateVecs(A, &x, NULL)); 60 61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 Create time stepping solver context 63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 64 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 65 PetscCall(TSSetType(ts, TSEIMEX)); 66 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &user)); 67 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user)); 68 PetscCall(TSSetIJacobian(ts, A, A, FormIJacobian, &user)); 69 ftime = 22; 70 PetscCall(TSSetMaxTime(ts, ftime)); 71 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 72 73 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 74 Set initial conditions 75 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 76 PetscCall(FormInitialSolution(ts, x, &user)); 77 PetscCall(TSSetSolution(ts, x)); 78 PetscCall(VecGetSize(x, &mx)); 79 80 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 81 Set runtime options 82 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 83 PetscCall(TSSetFromOptions(ts)); 84 85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86 Solve nonlinear system 87 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 88 PetscCall(TSSolve(ts, x)); 89 PetscCall(TSGetTime(ts, &ftime)); 90 PetscCall(TSGetStepNumber(ts, &steps)); 91 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "eps %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.param, steps, (double)ftime)); 92 /* PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));*/ 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Free work space. 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 PetscCall(MatDestroy(&A)); 98 PetscCall(VecDestroy(&x)); 99 PetscCall(TSDestroy(&ts)); 100 PetscCall(PetscFinalize()); 101 return 0; 102 } 103 104 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 105 AppCtx *user = (AppCtx *)ptr; 106 PetscScalar *f; 107 const PetscScalar *x; 108 PetscInt i, mx; 109 PetscReal hx, eps; 110 111 PetscFunctionBegin; 112 mx = user->mx; 113 eps = user->param; 114 hx = (user->xright - user->xleft) / (mx - 1); 115 PetscCall(VecGetArrayRead(X, &x)); 116 PetscCall(VecGetArray(F, &f)); 117 f[0] = 2. * eps * (x[1] - x[0]) / (hx * hx); /*boundary*/ 118 for (i = 1; i < mx - 1; i++) f[i] = eps * (x[i + 1] - 2. * x[i] + x[i - 1]) / (hx * hx); 119 f[mx - 1] = 2. * eps * (x[mx - 2] - x[mx - 1]) / (hx * hx); /*boundary*/ 120 PetscCall(VecRestoreArrayRead(X, &x)); 121 PetscCall(VecRestoreArray(F, &f)); 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr) { 126 AppCtx *user = (AppCtx *)ptr; 127 PetscScalar *f; 128 const PetscScalar *x, *xdot; 129 PetscInt i, mx; 130 131 PetscFunctionBegin; 132 mx = user->mx; 133 PetscCall(VecGetArrayRead(X, &x)); 134 PetscCall(VecGetArrayRead(Xdot, &xdot)); 135 PetscCall(VecGetArray(F, &f)); 136 137 for (i = 0; i < mx; i++) f[i] = xdot[i] - x[i] * (1. - x[i] * x[i]); 138 139 PetscCall(VecRestoreArrayRead(X, &x)); 140 PetscCall(VecRestoreArrayRead(Xdot, &xdot)); 141 PetscCall(VecRestoreArray(F, &f)); 142 PetscFunctionReturn(0); 143 } 144 145 static PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *ctx) { 146 AppCtx *user = (AppCtx *)ctx; 147 PetscScalar v; 148 const PetscScalar *x; 149 PetscInt i, col; 150 151 PetscFunctionBegin; 152 PetscCall(VecGetArrayRead(U, &x)); 153 for (i = 0; i < user->mx; i++) { 154 v = a - 1. + 3. * x[i] * x[i]; 155 col = i; 156 PetscCall(MatSetValues(J, 1, &i, 1, &col, &v, INSERT_VALUES)); 157 } 158 PetscCall(VecRestoreArrayRead(U, &x)); 159 160 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 161 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 162 if (J != Jpre) { 163 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 164 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 165 } 166 /* MatView(J,PETSC_VIEWER_STDOUT_WORLD);*/ 167 PetscFunctionReturn(0); 168 } 169 170 static PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ctx) { 171 AppCtx *user = (AppCtx *)ctx; 172 PetscInt i; 173 PetscScalar *x; 174 PetscReal hx, x_map; 175 176 PetscFunctionBegin; 177 hx = (user->xright - user->xleft) / (PetscReal)(user->mx - 1); 178 PetscCall(VecGetArray(U, &x)); 179 for (i = 0; i < user->mx; i++) { 180 x_map = user->xleft + i * hx; 181 if (x_map >= 0.7065) { 182 x[i] = PetscTanhReal((x_map - 0.8) / (2. * PetscSqrtReal(user->param))); 183 } else if (x_map >= 0.4865) { 184 x[i] = PetscTanhReal((0.613 - x_map) / (2. * PetscSqrtReal(user->param))); 185 } else if (x_map >= 0.28) { 186 x[i] = PetscTanhReal((x_map - 0.36) / (2. * PetscSqrtReal(user->param))); 187 } else if (x_map >= -0.7) { 188 x[i] = PetscTanhReal((0.2 - x_map) / (2. * PetscSqrtReal(user->param))); 189 } else if (x_map >= -1) { 190 x[i] = PetscTanhReal((x_map + 0.9) / (2. * PetscSqrtReal(user->param))); 191 } 192 } 193 PetscCall(VecRestoreArray(U, &x)); 194 PetscFunctionReturn(0); 195 } 196 197 /*TEST 198 199 test: 200 args: -ts_rtol 1e-04 -ts_dt 0.025 -pc_type lu -ksp_error_if_not_converged TRUE -ts_type eimex -ts_adapt_type none -ts_eimex_order_adapt -ts_eimex_max_rows 7 -ts_monitor_draw_solution 201 requires: x 202 timeoutfactor: 3 203 204 TEST*/ 205