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