|
| 1 | +## Solving Nonlinear PDEs with Devito {#sec-nonlin-devito} |
| 2 | + |
| 3 | +Having established the finite difference discretization of nonlinear PDEs, |
| 4 | +we now implement several solvers using Devito. The symbolic approach allows |
| 5 | +us to express nonlinear equations and handle the time-lagged coefficients |
| 6 | +naturally. |
| 7 | + |
| 8 | +### Nonlinear Diffusion: The Explicit Scheme |
| 9 | + |
| 10 | +The nonlinear diffusion equation |
| 11 | +$$ |
| 12 | +u_t = \nabla \cdot (D(u) \nabla u) |
| 13 | +$$ |
| 14 | +with solution-dependent diffusivity $D(u)$ requires special treatment. |
| 15 | +In 1D, the equation becomes: |
| 16 | +$$ |
| 17 | +u_t = \frac{\partial}{\partial x}\left(D(u) \frac{\partial u}{\partial x}\right) |
| 18 | +$$ |
| 19 | + |
| 20 | +For explicit time stepping, we evaluate $D$ at the previous time level: |
| 21 | +$$ |
| 22 | +u^{n+1}_i = u^n_i + \frac{\Delta t}{\Delta x^2} |
| 23 | +\left[D^n_{i+1/2}(u^n_{i+1} - u^n_i) - D^n_{i-1/2}(u^n_i - u^n_{i-1})\right] |
| 24 | +$$ |
| 25 | + |
| 26 | +where $D^n_{i+1/2} = \frac{1}{2}(D(u^n_i) + D(u^n_{i+1}))$. |
| 27 | + |
| 28 | +### The Devito Implementation |
| 29 | + |
| 30 | +```python |
| 31 | +from devito import Grid, TimeFunction, Eq, Operator, Constant |
| 32 | +import numpy as np |
| 33 | + |
| 34 | +# Domain and discretization |
| 35 | +L = 1.0 # Domain length |
| 36 | +Nx = 100 # Grid points |
| 37 | +T = 0.1 # Final time |
| 38 | +F = 0.4 # Target Fourier number |
| 39 | + |
| 40 | +dx = L / Nx |
| 41 | +D_max = 1.0 # Maximum diffusion coefficient |
| 42 | +dt = F * dx**2 / D_max # Time step from stability |
| 43 | + |
| 44 | +# Create Devito grid |
| 45 | +grid = Grid(shape=(Nx + 1,), extent=(L,)) |
| 46 | + |
| 47 | +# Time-varying field with space_order=2 for halo access |
| 48 | +u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2) |
| 49 | +``` |
| 50 | + |
| 51 | +### Handling the Nonlinear Diffusion Coefficient |
| 52 | + |
| 53 | +For nonlinear diffusion, the diffusivity depends on the solution. Common |
| 54 | +forms include: |
| 55 | + |
| 56 | +| Type | $D(u)$ | Application | |
| 57 | +|------|--------|-------------| |
| 58 | +| Constant | $D_0$ | Linear heat conduction | |
| 59 | +| Linear | $D_0(1 + \alpha u)$ | Temperature-dependent conductivity | |
| 60 | +| Porous medium | $D_0 m u^{m-1}$ | Flow in porous media | |
| 61 | + |
| 62 | +The `src.nonlin` module provides several diffusion coefficient functions: |
| 63 | + |
| 64 | +```python |
| 65 | +from src.nonlin import ( |
| 66 | + constant_diffusion, |
| 67 | + linear_diffusion, |
| 68 | + porous_medium_diffusion, |
| 69 | +) |
| 70 | + |
| 71 | +# Constant D(u) = 1.0 |
| 72 | +D_const = lambda u: constant_diffusion(u, D0=1.0) |
| 73 | + |
| 74 | +# Linear D(u) = 1 + 0.5*u |
| 75 | +D_linear = lambda u: linear_diffusion(u, D0=1.0, alpha=0.5) |
| 76 | + |
| 77 | +# Porous medium D(u) = 2*u (m=2) |
| 78 | +D_porous = lambda u: porous_medium_diffusion(u, m=2.0, D0=1.0) |
| 79 | +``` |
| 80 | + |
| 81 | +### Complete Nonlinear Diffusion Solver |
| 82 | + |
| 83 | +The `src.nonlin` module provides `solve_nonlinear_diffusion_explicit`: |
| 84 | + |
| 85 | +```python |
| 86 | +from src.nonlin import solve_nonlinear_diffusion_explicit |
| 87 | +import numpy as np |
| 88 | + |
| 89 | +# Initial condition: smooth bump |
| 90 | +def I(x): |
| 91 | + return np.sin(np.pi * x) |
| 92 | + |
| 93 | +result = solve_nonlinear_diffusion_explicit( |
| 94 | + L=1.0, # Domain length |
| 95 | + Nx=100, # Grid points |
| 96 | + T=0.1, # Final time |
| 97 | + F=0.4, # Fourier number |
| 98 | + I=I, # Initial condition |
| 99 | + D_func=lambda u: linear_diffusion(u, D0=1.0, alpha=0.5), |
| 100 | +) |
| 101 | + |
| 102 | +print(f"Final time: {result.t:.4f}") |
| 103 | +print(f"Max solution: {result.u.max():.6f}") |
| 104 | +``` |
| 105 | + |
| 106 | +### Reaction-Diffusion with Operator Splitting |
| 107 | + |
| 108 | +The reaction-diffusion equation |
| 109 | +$$ |
| 110 | +u_t = a u_{xx} + R(u) |
| 111 | +$$ |
| 112 | +combines diffusion with a nonlinear reaction term. Operator splitting |
| 113 | +separates these effects: |
| 114 | + |
| 115 | +**Lie Splitting (first-order):** |
| 116 | +1. Solve $u_t = a u_{xx}$ for time $\Delta t$ |
| 117 | +2. Solve $u_t = R(u)$ for time $\Delta t$ |
| 118 | + |
| 119 | +**Strang Splitting (second-order):** |
| 120 | +1. Solve $u_t = R(u)$ for time $\Delta t/2$ |
| 121 | +2. Solve $u_t = a u_{xx}$ for time $\Delta t$ |
| 122 | +3. Solve $u_t = R(u)$ for time $\Delta t/2$ |
| 123 | + |
| 124 | +### Reaction Terms |
| 125 | + |
| 126 | +The module provides common reaction terms: |
| 127 | + |
| 128 | +```python |
| 129 | +from src.nonlin import ( |
| 130 | + logistic_reaction, |
| 131 | + fisher_reaction, |
| 132 | + allen_cahn_reaction, |
| 133 | +) |
| 134 | + |
| 135 | +# Logistic growth: R(u) = r*u*(1 - u/K) |
| 136 | +R_logistic = lambda u: logistic_reaction(u, r=1.0, K=1.0) |
| 137 | + |
| 138 | +# Fisher-KPP: R(u) = r*u*(1 - u) |
| 139 | +R_fisher = lambda u: fisher_reaction(u, r=1.0) |
| 140 | + |
| 141 | +# Allen-Cahn: R(u) = u - u^3 |
| 142 | +R_allen_cahn = lambda u: allen_cahn_reaction(u, epsilon=1.0) |
| 143 | +``` |
| 144 | + |
| 145 | +### Reaction-Diffusion Solver |
| 146 | + |
| 147 | +```python |
| 148 | +from src.nonlin import solve_reaction_diffusion_splitting |
| 149 | + |
| 150 | +# Initial condition with small perturbation |
| 151 | +def I(x): |
| 152 | + return 0.5 * np.sin(np.pi * x) |
| 153 | + |
| 154 | +# Strang splitting (second-order) |
| 155 | +result = solve_reaction_diffusion_splitting( |
| 156 | + L=1.0, |
| 157 | + a=0.1, # Diffusion coefficient |
| 158 | + Nx=100, |
| 159 | + T=0.5, |
| 160 | + F=0.4, |
| 161 | + I=I, |
| 162 | + R_func=lambda u: fisher_reaction(u, r=1.0), |
| 163 | + splitting="strang", |
| 164 | +) |
| 165 | +``` |
| 166 | + |
| 167 | +The Strang splitting achieves second-order accuracy in time, while Lie |
| 168 | +splitting is only first-order. For problems with fast reactions or |
| 169 | +long simulation times, the higher accuracy of Strang splitting is |
| 170 | +beneficial. |
| 171 | + |
| 172 | +### Burgers' Equation |
| 173 | + |
| 174 | +The viscous Burgers' equation |
| 175 | +$$ |
| 176 | +u_t + u u_x = \nu u_{xx} |
| 177 | +$$ |
| 178 | +is a prototype for nonlinear advection with viscous dissipation. The |
| 179 | +nonlinear term $u u_x$ can cause shock formation for small $\nu$. |
| 180 | + |
| 181 | +We use the conservative form $(u^2/2)_x$ with centered differences: |
| 182 | + |
| 183 | +```python |
| 184 | +from src.nonlin import solve_burgers_equation |
| 185 | + |
| 186 | +result = solve_burgers_equation( |
| 187 | + L=2.0, # Domain length |
| 188 | + nu=0.01, # Viscosity |
| 189 | + Nx=100, # Grid points |
| 190 | + T=0.5, # Final time |
| 191 | + C=0.5, # Target CFL number |
| 192 | +) |
| 193 | +``` |
| 194 | + |
| 195 | +### Stability for Burgers' Equation |
| 196 | + |
| 197 | +The time step must satisfy both the CFL condition for advection: |
| 198 | +$$ |
| 199 | +C = \frac{|u|_{\max} \Delta t}{\Delta x} \le 1 |
| 200 | +$$ |
| 201 | + |
| 202 | +and the diffusion stability condition: |
| 203 | +$$ |
| 204 | +F = \frac{\nu \Delta t}{\Delta x^2} \le 0.5 |
| 205 | +$$ |
| 206 | + |
| 207 | +The solver automatically chooses $\Delta t$ to satisfy both conditions |
| 208 | +with a safety factor. |
| 209 | + |
| 210 | +### The Effect of Viscosity |
| 211 | + |
| 212 | +```python |
| 213 | +import matplotlib.pyplot as plt |
| 214 | + |
| 215 | +fig, axes = plt.subplots(1, 2, figsize=(12, 5)) |
| 216 | + |
| 217 | +for ax, nu in zip(axes, [0.1, 0.01]): |
| 218 | + result = solve_burgers_equation( |
| 219 | + L=2.0, nu=nu, Nx=100, T=0.5, C=0.3, |
| 220 | + I=lambda x: np.sin(np.pi * x), |
| 221 | + save_history=True, |
| 222 | + ) |
| 223 | + |
| 224 | + for i in range(0, len(result.t_history), len(result.t_history)//5): |
| 225 | + ax.plot(result.x, result.u_history[i], |
| 226 | + label=f't = {result.t_history[i]:.2f}') |
| 227 | + |
| 228 | + ax.set_xlabel('x') |
| 229 | + ax.set_ylabel('u') |
| 230 | + ax.set_title(f'Burgers, nu = {nu}') |
| 231 | + ax.legend() |
| 232 | +``` |
| 233 | + |
| 234 | +Higher viscosity ($\nu = 0.1$) smooths the solution, while lower |
| 235 | +viscosity ($\nu = 0.01$) allows steeper gradients to develop. |
| 236 | + |
| 237 | +### Picard Iteration for Implicit Schemes |
| 238 | + |
| 239 | +For stiff nonlinear problems, implicit time stepping may be necessary. |
| 240 | +Picard iteration solves the nonlinear system by repeated linearization: |
| 241 | + |
| 242 | +1. Guess $u^{n+1, (0)} = u^n$ |
| 243 | +2. For $k = 0, 1, 2, \ldots$: |
| 244 | + - Evaluate $D^{(k)} = D(u^{n+1, (k)})$ |
| 245 | + - Solve the linear system for $u^{n+1, (k+1)}$ |
| 246 | + - Check convergence: $\|u^{n+1, (k+1)} - u^{n+1, (k)}\| < \epsilon$ |
| 247 | + |
| 248 | +```python |
| 249 | +from src.nonlin import solve_nonlinear_diffusion_picard |
| 250 | + |
| 251 | +result = solve_nonlinear_diffusion_picard( |
| 252 | + L=1.0, |
| 253 | + Nx=50, |
| 254 | + T=0.05, |
| 255 | + dt=0.001, # Can use larger dt than explicit |
| 256 | +) |
| 257 | +``` |
| 258 | + |
| 259 | +The implicit scheme removes the time step restriction but requires |
| 260 | +solving a linear system at each iteration. |
| 261 | + |
| 262 | +### Summary |
| 263 | + |
| 264 | +Key points for nonlinear PDEs with Devito: |
| 265 | + |
| 266 | +1. **Nonlinear diffusion**: Use explicit scheme with lagged coefficient |
| 267 | + evaluation and Fourier number $F \le 0.5$ |
| 268 | +2. **Operator splitting**: Separates diffusion and reaction for |
| 269 | + reaction-diffusion equations; Strang is second-order |
| 270 | +3. **Burgers' equation**: Requires both CFL and diffusion stability |
| 271 | + conditions; viscosity controls smoothness |
| 272 | +4. **Picard iteration**: Enables implicit schemes for stiff problems |
| 273 | + at the cost of solving linear systems |
| 274 | + |
| 275 | +The `src.nonlin` module provides: |
| 276 | +- `solve_nonlinear_diffusion_explicit` |
| 277 | +- `solve_reaction_diffusion_splitting` |
| 278 | +- `solve_burgers_equation` |
| 279 | +- `solve_nonlinear_diffusion_picard` |
| 280 | +- Diffusion coefficients: `constant_diffusion`, `linear_diffusion`, |
| 281 | + `porous_medium_diffusion` |
| 282 | +- Reaction terms: `logistic_reaction`, `fisher_reaction`, |
| 283 | + `allen_cahn_reaction` |
0 commit comments