Add diffusion-implicit heat equation tutorial #476
Conversation
| Next, we'll define some global problem parameters: | ||
| ```julia | ||
| a,b, n = 0, 1, 10 # zmin, zmax, number of cells | ||
| n̂_min, n̂_max = -1, 1 # Outward facing unit vectors | ||
| α = 100; # thermal diffusivity, larger means more stiff | ||
| β, γ = 10000, π; # source term coefficients | ||
| Δt = 1000; # timestep size | ||
| N_t = 10; # number of timesteps to take | ||
| FT = Float64; # float type | ||
| Δz = FT(b-a)/FT(n) | ||
| Δz² = Δz^2; | ||
| ∇²_op = [1/Δz², -2/Δz², 1/Δz²]; # interior Laplacian operator | ||
| ∇T_bottom = 10; # Temperature gradient at the top | ||
| T_top = 1; # Temperature at the bottom | ||
| S(z) = β*sin(γ*z) # source term, (sin for easy integration) | ||
| zf = range(a, b, length=n+1); # coordinates on cell faces |
There was a problem hiding this comment.
Is this actually needed for the tutorial? It seems like this is doing a lot of unnecessary things. Seems more like a https://github.com/SciML/SciMLTutorials.jl tutorial instead of something streamlined for the docs.
There was a problem hiding this comment.
Ah, yeah, I wasn’t exactly sure where to put this. I’m fine with moving it (the whole tutorial, right?) there
There was a problem hiding this comment.
Yeah, let's go with the move. I am putting up a bounty soon anyways to try and find a way to incorporate those tutorials and the benchmarks into the documentation. It's mostly split so that it's easy to regenerate the results using the AMDCI computers, allowing you to go ham and make it a big equation and really demonstrate sparsity handling etc. I hope that it will soon then feed back nicely into the main docs.
| @@ -0,0 +1,186 @@ | |||
| # Solving the heat equation with diffusion-implicit time-stepping. | |||
There was a problem hiding this comment.
| # Solving the heat equation with diffusion-implicit time-stepping. | |
| # IMEX Integrators |
This title is more informative as to what the tutorial is actually trying to teach.
There was a problem hiding this comment.
If it's moved, the original title makes sense though.
| ```math | ||
| \frac{∂²T}{∂²z} = -S(z)/α = -β \sin(γ z)/α \\ | ||
| \frac{∂T}{∂z} = β \cos(γ z)/(γ α)+c_1 \\ | ||
| T(z) = β \sin(γ z)/(γ^2 α)+c_1 z+c_2, \qquad \text{(generic solution)} | ||
| ``` | ||
| Apply bottom boundary condition: | ||
| ```math | ||
| \frac{∂T}{∂z}(a) = β \cos(γ a)/(γ α)+c_1 = ∇T_{bottom} \\ | ||
| c_1 = ∇T_{bottom}-β \cos(γ a)/(γ α) | ||
| ``` | ||
|
|
||
| Apply top boundary condition: | ||
| ```math | ||
| T(b) = β \sin(γ b)/(γ^2 α)+c_1 b+c_2 = T_{top} \\ | ||
| c_2 = T_{top}-(β \sin(γ b)/(γ^2 α)+c_1 b) | ||
| ``` | ||
|
|
||
| And now let's define this in a julia function: | ||
| ```julia | ||
| function T_analytic(z) # Analytic steady state solution | ||
| c1 = ∇T_bottom-β*cos(γ*a)/(γ*α) | ||
| c2 = T_top-(β*sin(γ*b)/(γ^2*α)+c1*b) | ||
| return β*sin(γ*z)/(γ^2*α)+c1*z+c2 | ||
| end | ||
| ``` | ||
|
|
||
| Here, we'll derivation the matrix form of the temporal discretization we wish to use (diffusion-implicit and explicit Euler): | ||
| ```math | ||
| ∂_t T = α ∇²T + S \\ | ||
| (T^{n+1}-T^n) = Δt (α ∇²T^{n+1} + S) \\ | ||
| (T^{n+1} - Δt α ∇²T^{n+1}) = T^n + Δt S \\ | ||
| (I - Δt α ∇²) T^{n+1} = T^n + Δt S | ||
| ``` | ||
|
|
||
| Using the following index denotion: | ||
|
|
||
| - `i` first interior index | ||
| - `b` boundary index | ||
| - `g` ghost index | ||
|
|
||
| we'll derive the Dirichlet boundary stencil & source: | ||
| ```math | ||
| ∂_t T = α (T[i-1]+T[b]-2 T[i])/Δz² + S \\ | ||
| ∂_t T = α (T[i-1]-2 T[i])/Δz² + S + α T[b] / Δz² | ||
| ``` | ||
|
|
||
| and Neumann boundary stencil & source: | ||
| ```math | ||
| ∇T_{bottom} n̂ = (T[g] - T[i])/(2Δz), \qquad n̂ = [-1,1] ∈ [zmin,zmax] \\ | ||
| T[i] + 2 Δz ∇T_{bottom} n̂ = T[g] \\ | ||
| ∂_t T = α (((T[i] + 2 Δz ∇T_{bottom} n̂) - T[b])/Δz - (T[b] - T[i])/Δz)/Δz + S \\ | ||
| ∂_t T = α (((T[i]) - T[b])/Δz - (T[b] - T[i])/Δz)/Δz + S + α 2 Δz ∇T_{bottom} n̂/Δz² \\ | ||
| ∂_t T = α (2 T[i] - 2 T[b])/Δz² + S + 2α ∇T_{bottom} n̂/Δz | ||
| ``` |
There was a problem hiding this comment.
Yeah, all of this is too much build up and loses the message. I think this should all go to https://github.com/SciML/SciMLTutorials.jl . It would be good to get an IMEX in here, but such a tutorial would want to have no more than a 1 paragraph problem statement, which is very different from what this is trying to do. This is trying to show a workflow specific to a use cases, which is what the extended tutorial set is made for.
| tspan, | ||
| params | ||
| ) | ||
| alg = IMEXEuler(linsolve=LinSolveFactorize(lu!)) |
There was a problem hiding this comment.
Are you sure IMEXEuler is a good idea here? I would assume KenCarp4 etc. would perform better.
There was a problem hiding this comment.
Yeah, probably, I was just trying to keep certain things simple.
|
Closing in favor of SciMLTutorials#416 |
2 participants
This PR adds a diffusion-implicit heat tutorial. cc @ChrisRackauckas