2D with parallel map
For this example we will solve the head equation with a mangetic field aligned with the grid
\[ \mathbf{B} = (0,0,1)\]
In this case we expect the parallel operator to do nothing since $\mathbf{P}_f=\mathbf{P}_b=I$
using FaADEFor this we'll solve the field aligned equation is
\[ \frac{\partial u}{\partial t} = \kappa_\perp \nabla_\perp^2 u + \mathcal{P}_\parallel u\]
with Dirichlet boundaries in $x$
\[ u(0,y,t) = 0, \qquad \partial_x u(1,y,t) = 0\]
periodic in $y$, and initial condition
\[ u(x,0) = \exp\left(\frac{-(x-0.5)^2}{0.02}\right)\]
We first need to create a domain to solve the PDE using Grid2D
𝒟x = [0.0,1.0]
𝒟y = [0.0,1.0]
nx = ny = 21
grid = Grid2D(𝒟x,𝒟y,nx,ny)Grid2D{Float64, CartesianMetric, Matrix{Float64}}([0.0 0.0 … 0.0 0.0; 0.05 0.05 … 0.05 0.05; … ; 0.95 0.95 … 0.95 0.95; 1.0 1.0 … 1.0 1.0], [0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0; … ; 0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0], 0.05, 0.05, 21, 21, [1.0;;], [1.0;;], [1.0;;], [1.0;;], [1.0;;])The initial condition is
u₀(x,y) = exp(-((x-0.5)^2 + (y-0.5)^2) / 0.02)u₀ (generic function with 1 method)Set the FD order to 2,
order = 22The boundary conditions are defined by creating SimultanousApproximationTerm objects, which will then be fed to the PDE structure
BoundaryLeft = SAT_Dirichlet((X,t)->0.0,grid.Δx, Left, order)
BoundaryRight= SAT_Neumann((X,t)->0.0, grid.Δx, Right, order)
BoundaryUp = SAT_Periodic(grid.Δy, Up, order)
BoundaryDown = SAT_Periodic(grid.Δy, Down, order)
BCs = (BoundaryLeft,BoundaryRight,BoundaryUp,BoundaryDown)(SAT_Dirichlet{FaADE.Helpers.NodeType{:Left, 1}, :Cartesian, Float64, Vector{Float64}, Main.var"#1#2", typeof(eachcol)}(FaADE.Helpers.NodeType{:Left, 1}(), 1, 2, Main.var"#1#2"(), 1600.0, [800.0, -400.0], 0.05, 1.0, [-2.0], eachcol, 0.0, :Cartesian), SAT_Neumann{FaADE.Helpers.NodeType{:Right, 1}, :Cartesian, Float64, Vector{Float64}, Main.var"#3#4", typeof(eachcol)}(FaADE.Helpers.NodeType{:Right, 1}(), 1, 2, Main.var"#3#4"(), -40.0, [-20.0, 20.0], 0.05, 1.0, eachcol, 0.0, :Cartesian), SAT_Periodic{FaADE.Helpers.NodeType{:Right, 2}, :Cartesian, Float64, Vector{Float64}, FaADE.SATs.var"#τ₀#1"{Float64, Float64}, typeof(eachrow)}(FaADE.Helpers.BoundaryConditionType{:Periodic}(), FaADE.Helpers.NodeType{:Right, 2}(), 2, 2, [-40.0, 20.0], [-20.0, 40.0], [-20.0, 20.0], [-20.0, 20.0], 0.05, 20.0, -10.0, FaADE.SATs.var"#τ₀#1"{Float64, Float64}(0.05, 0.5), eachrow, 0.0, :Cartesian), SAT_Periodic{FaADE.Helpers.NodeType{:Left, 2}, :Cartesian, Float64, Vector{Float64}, FaADE.SATs.var"#τ₀#1"{Float64, Float64}, typeof(eachrow)}(FaADE.Helpers.BoundaryConditionType{:Periodic}(), FaADE.Helpers.NodeType{:Left, 2}(), 2, 2, [-40.0, 20.0], [-20.0, 40.0], [-20.0, 20.0], [-20.0, 20.0], 0.05, 20.0, -10.0, FaADE.SATs.var"#τ₀#1"{Float64, Float64}(0.05, 0.5), eachrow, 0.0, :Cartesian))The 2 input to the periodic boundary ensures it is along the y-axis.
Set the diffusion in $x$ and $y$ directions to 1
Kx = Ky = 1.01.0The parallel penalty function can be generated by providing the code the ODE for magnetic field lines. In this case we'll just as a magnetic field that does nothing,
function Bfield(X,x,p,t)
X[2] = 0.0
X[1] = 0.0
endBfield (generic function with 1 method)Assuming a $2\pi$ periodicity then we can construct a parallel grid object with construct_grid, then generate the ParallelData for passing to the problem object,
PGrid = construct_grid(Bfield,grid,[-2π,2π])
PData = ParallelData(PGrid,grid,order,κ=1.0)ParallelData{Float64, 2, FaADE.ParallelOperator.ParallelGrid{Float64, 2, FaADE.ParallelOperator.ParGrid{Float64, Matrix{Float64}}, Matrix{Float64}}, FaADE.ParallelOperator.MagneticField{Nothing, :EQUILIBRIUM, Float64, Matrix{Float64}}, CubicHermiteSpline.BivariateCHSInterpolation{Float64, DelaunayTriangulation.Triangulation{Matrix{Float64}, Set{Tuple{Int64, Int64, Int64}}, Vector{Int64}, DelaunayTriangulation.ZeroWeight{Float64}, Int64, Tuple{Int64, Int64}, Set{Tuple{Int64, Int64}}, Tuple{}, Dict{Tuple{Int64, Int64}, Tuple{Vector{Int64}, Int64}}, Dict{Int64, Vector{Int64}}, Dict{Int64, UnitRange{Int64}}, Dict{Int64, DelaunayTriangulation.RepresentativeCoordinates{Int64, Float64}}, DelaunayTriangulation.TriangulationCache{DelaunayTriangulation.Triangulation{Matrix{Float64}, Set{Tuple{Int64, Int64, Int64}}, Vector{Int64}, DelaunayTriangulation.ZeroWeight{Float64}, Int64, Tuple{Int64, Int64}, Set{Tuple{Int64, Int64}}, Tuple{}, Dict{Tuple{Int64, Int64}, Tuple{Vector{Int64}, Int64}}, Dict{Int64, Vector{Int64}}, Dict{Int64, UnitRange{Int64}}, Dict{Int64, DelaunayTriangulation.RepresentativeCoordinates{Int64, Float64}}, DelaunayTriangulation.TriangulationCache{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Nothing}, Vector{Int64}, Set{Tuple{Int64, Int64}}, Vector{Int64}, Set{Tuple{Int64, Int64, Int64}}, NTuple{4, SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}}, Tuple{SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}, SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}}, NTuple{27, SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}}}, Nothing}}, Nothing}(FaADE.ParallelOperator.ParallelGrid{Float64, 2, FaADE.ParallelOperator.ParGrid{Float64, Matrix{Float64}}, Matrix{Float64}}(FaADE.ParallelOperator.ParGrid{Float64, Matrix{Float64}}([0.0 0.0 … 0.0 0.0; 0.05 0.05 … 0.05 0.05; … ; 0.95 0.95 … 0.95 0.95; 1.0 1.0 … 1.0 1.0], [0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0; … ; 0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0], [0;;]), FaADE.ParallelOperator.ParGrid{Float64, Matrix{Float64}}([0.0 0.0 … 0.0 0.0; 0.05 0.05 … 0.05 0.05; … ; 0.95 0.95 … 0.95 0.95; 1.0 1.0 … 1.0 1.0], [0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0; … ; 0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0], [0;;])), 1.0, 1.0, nothing, Bivariate Cubic Hermite Spline Interpolation
Number of points: 441
Number of Delaunay triangles: 800
x: [0.0, 1.0], y: [0.0, 1.0]
f: [0.0, 0.0]
fx: [0.0, 0.0], fy: [0.0, 0.0]
, 0.05, 0.05, FaADE.Derivatives.CompositeH{2, Float64, Vector{Float64}, FaADE.Derivatives.DiagonalH{Float64, Vector{Float64}}}((FaADE.Derivatives.DiagonalH{Float64, Vector{Float64}}([0.5], 1.0, 0.05, 21, 1), FaADE.Derivatives.DiagonalH{Float64, Vector{Float64}}([0.5], 1.0, 0.05, 21, 1)), [21, 21]), [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0;;], FaADE.ParallelOperator.MagneticField{Nothing, :EQUILIBRIUM, Float64, Matrix{Float64}}(nothing, false, [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]), [0.0])Lastly before solving we define our time step and simulation time,
Δt = 0.01grid.Δx;
t_f = 100Δt;Now we can create a Problem2D object to pass to the solver,
P = Problem2D(order,u₀,Kx,Ky,grid,BCs,parallel=PData)2 dimensional PDE ProblemFinally we call the solver
soln = solve(P,grid,Δt,t_f)FaADE.solvers.solution{Float64, Matrix{Float64}, Grid2D{Float64, CartesianMetric, Matrix{Float64}}, Problem2D{Float64, 2, Float64, FaADE.Helpers.SourceTerm{Nothing}, Tuple{SAT_Dirichlet{FaADE.Helpers.NodeType{:Left, 1}, :Cartesian, Float64, Vector{Float64}, Main.var"#1#2", typeof(eachcol)}, SAT_Neumann{FaADE.Helpers.NodeType{:Right, 1}, :Cartesian, Float64, Vector{Float64}, Main.var"#3#4", typeof(eachcol)}, SAT_Periodic{FaADE.Helpers.NodeType{:Right, 2}, :Cartesian, Float64, Vector{Float64}, FaADE.SATs.var"#τ₀#1"{Float64, Float64}, typeof(eachrow)}, SAT_Periodic{FaADE.Helpers.NodeType{:Left, 2}, :Cartesian, Float64, Vector{Float64}, FaADE.SATs.var"#τ₀#1"{Float64, Float64}, typeof(eachrow)}}, ParallelData{Float64, 2, FaADE.ParallelOperator.ParallelGrid{Float64, 2, FaADE.ParallelOperator.ParGrid{Float64, Matrix{Float64}}, Matrix{Float64}}, FaADE.ParallelOperator.MagneticField{Nothing, :EQUILIBRIUM, Float64, Matrix{Float64}}, CubicHermiteSpline.BivariateCHSInterpolation{Float64, DelaunayTriangulation.Triangulation{Matrix{Float64}, Set{Tuple{Int64, Int64, Int64}}, Vector{Int64}, DelaunayTriangulation.ZeroWeight{Float64}, Int64, Tuple{Int64, Int64}, Set{Tuple{Int64, Int64}}, Tuple{}, Dict{Tuple{Int64, Int64}, Tuple{Vector{Int64}, Int64}}, Dict{Int64, Vector{Int64}}, Dict{Int64, UnitRange{Int64}}, Dict{Int64, DelaunayTriangulation.RepresentativeCoordinates{Int64, Float64}}, DelaunayTriangulation.TriangulationCache{DelaunayTriangulation.Triangulation{Matrix{Float64}, Set{Tuple{Int64, Int64, Int64}}, Vector{Int64}, DelaunayTriangulation.ZeroWeight{Float64}, Int64, Tuple{Int64, Int64}, Set{Tuple{Int64, Int64}}, Tuple{}, Dict{Tuple{Int64, Int64}, Tuple{Vector{Int64}, Int64}}, Dict{Int64, Vector{Int64}}, Dict{Int64, UnitRange{Int64}}, Dict{Int64, DelaunayTriangulation.RepresentativeCoordinates{Int64, Float64}}, DelaunayTriangulation.TriangulationCache{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Nothing}, Vector{Int64}, Set{Tuple{Int64, Int64}}, Vector{Int64}, Set{Tuple{Int64, Int64, Int64}}, NTuple{4, SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}}, Tuple{SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}, SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}}, NTuple{27, SubArray{Float64, 1, Memory{Float64}, Tuple{UnitRange{Int64}}, true}}}, Nothing}}, Nothing}}}([[1.3887943864964021e-11 1.493094676197164e-10 … 1.4930946761971692e-10 1.3887943864964021e-11; 1.493094676197164e-10 1.6052280551856116e-9 … 1.6052280551856172e-9 1.493094676197164e-10; … ; 1.4930946761971692e-10 1.6052280551856172e-9 … 1.6052280551856172e-9 1.4930946761971692e-10; 1.3887943864964021e-11 1.493094676197164e-10 … 1.4930946761971692e-10 1.3887943864964021e-11], [2.188265059742769e-27 -2.1853305122948477e-27 … -2.185330512289783e-27 2.188265059742769e-27; 0.008437955927116107 0.008561560976077987 … 0.008561560976077985 0.008437955927116107; … ; 0.03765454097133522 0.03820620389495539 … 0.038206203894955396 0.03765454097133522; 0.037112829058002256 0.03765655996599567 … 0.03765655996599567 0.037112829058002256]], Grid2D{Float64, CartesianMetric, Matrix{Float64}}([0.0 0.0 … 0.0 0.0; 0.05 0.05 … 0.05 0.05; … ; 0.95 0.95 … 0.95 0.95; 1.0 1.0 … 1.0 1.0], [0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0; … ; 0.0 0.05 … 0.95 1.0; 0.0 0.05 … 0.95 1.0], 0.05, 0.05, 21, 21, [1.0;;], [1.0;;], [1.0;;], [1.0;;], [1.0;;]), [0.0005, 0.0005], [0.0, 0.05], 2 dimensional PDE Problem, 0.0053009123599797574, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])We can now plot the solution
using Plots
surface(grid.gridx, grid.gridy ,soln.u[2])
This page was generated using Literate.jl.