2D Julia Operators

julia/MOLE.jl/docs/make.jl

@@ -1,5 +1,5 @@
 push!(LOAD_PATH,"../src/")
-using Documenter, MOLE
+using Documenter, MOLE, SparseArrays
 
 makedocs(
     sitename="MOLE.jl Docs",

julia/MOLE.jl/docs/src/index.md

@@ -45,6 +45,7 @@ julia --project=. path/to/myScript.jl
 To run the unit tests, first enter the Julia REPL as in the above section (that is, by running the command `julia --project=.` from the directory `mole/julia/MOLE.jl`). Next, enter the `pkg` mode by pressing `]`, then type the command `test`. The results of the unit tests should be displayed to your console.
 
 ## Building the documentation
+
 MOLE.jl uses [Documenter.jl](https://documenter.juliadocs.org/stable/) to build its Julia implementation documentation. From the `mole/julia/MOLE.jl` directory, navigate to the `docs/` directory, with 
 
 ```sh
@@ -62,6 +63,7 @@ include("make.jl")
 ```
 
 ### From the command line
+
 To build the documentation from the `docs/` directory, use the following commands to instatiate and precompile the documentation environment:
 
 ```sh
@@ -79,22 +81,35 @@ julia --project=. make.jl
 ```
 
 ### Preview the documentation
+
 Once you have built the documentation (either from the REPL or the command line), you can inspect the documentation in `docs/build` with the `index.html` file as the homepage.
 
 ## Functions
 
 ### Operators
 
 ```@docs
-MOLE.div(k::Int,m::Int,dx)
-MOLE.grad(k::Int,m::Int,dx)
-MOLE.lap(k::Int,m::Int,dx)
+MOLE.Operators.div(k::Int, m::Int, dx; dc::NTuple{2,T}, nc::NTuple{2,T})
+MOLE.Operators.div(k::Int, ticks::AbstractVector)
+MOLE.Operators.div(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T}, nc::NTuple{4,T})
+MOLE.Operators.div(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+MOLE.Operators.grad(k::Int, m::Int, dx; dc::NTuple{2,T}, nc::NTuple{2,T})
+MOLE.Operators.grad(k::Int, ticks::AbstractVector)
+MOLE.Operators.grad(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T}, nc::NTuple{4,T})
+MOLE.Operators.grad(k::Int, xticks::AbstractVector, yticks::AbstractVector)
+MOLE.Operators.lap(k::Int, m::Int, dx; dc::NTuple{2,T}, nc::NTuple{2,T})
+MOLE.Operators.lap(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T}, nc::NTuple{4,T})
 ```
 
 ### Utilities
 
 ```@docs
-MOLE.robinBC(k::Int, m::Int, dx, a, b)
+MOLE.BCs.robinBC(k::Int, m::Int, dx, a, b)
+MOLE.BCs.robinBC(k::Int, m::Int, dx, n::Int, dy, a, b)
+MOLE.BCs.ScalarBC1D{T}
+MOLE.BCs.ScalarBC2D{T}
+MOLE.BCs.addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::MOLE.BCs.ScalarBC1D, k::Integer, m::Integer, dx)
+MOLE.BCs.addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::MOLE.BCs.ScalarBC2D, k::Integer, m::Integer, dx, n::Integer, dy)
 ```
 
 ## Examples
@@ -104,5 +119,11 @@ The MOLE library contains examples demonstrating how to use the operators, in a
 Currently, the following examples are available in the MOLE Julia package.
 
 - Elliptic Problems
-    - 1D Examples
-        - ```elliptic1D```: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators.
+  - 1D Examples
+    - ```elliptic1D```: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators.
+  - 2D Examples
+    - ```elliptic2DXDirichletYDirichlet```: A script that solves the 2D Laplace equation, $\nabla^2 u = 0$, with Dirichlet boundary conditions in $x$ and $y$ using mimetic operators.
+    - ```elliptic2DXPerYDirichlet```: A script that solves the 2D Laplace equation, $\nabla^2 u = 0$, with periodic bonudary conditions in $x$ and Dirichlet boundary conditions in $y$ using mimetic operators.
+- Parabolic Problems
+  - 2D Examples
+    - ```parabolic2D```: A script that solves the 2D heat equation, $u_t = \nu \nabla^2 u$, with Dirichlet boundary conditions in $x$ and $y$ using mimetic operators.

julia/MOLE.jl/examples/elliptic2DXDirichletYDirichlet.jl

@@ -0,0 +1,56 @@
+using LinearAlgebra
+using SparseArrays
+using Plots
+
+import MOLE: Operators, BCs
+
+# Parameters
+k = 2       # Order of accuracy
+m = 99      # Number of cells in x-direction
+n = m + 2   # Number of cells in y-direction
+dx = pi / m # Step size in x-direction
+dy = pi / n # Step size in y-direction
+
+# Grid
+xc = [0; (dx / 2.0) : dx : (pi - dx / 2.0); pi]
+yc = [0; (dy / 2.0) : dy : (pi - dy / 2.0); pi]
+X = ones(1, n + 2) .* xc
+Y = yc' .* ones(m + 2, 1)
+
+# Exact Solution
+ue = exp.(X) .* cos.(Y)
+
+# Boundary Conditions
+dc = (1.0, 1.0, 1.0, 1.0)
+nc = (0.0, 0.0, 0.0, 0.0)
+v = (vec(ue[1,2:end-1]'), vec(ue[end,2:end-1]'), vec(ue[:,1]), vec(ue[:, end]))
+bc = BCs.ScalarBC2D(dc, nc, v)
+
+# Operator
+A = - Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+
+# RHS
+b = zeros(m + 2, n + 2)
+b = vec(b)
+
+# Apply BCs
+A0, b0 = BCs.addScalarBC!(sparse(A), b, bc, k, m, dx, n, dy)
+
+# Approximate Solution
+ua = A0 \ b0
+ua = Matrix(reshape(ua, m + 2, n + 2))
+
+pa = heatmap(xc, yc, ua, title = "Approximate Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pa)
+println("Press Enter to close the plot and open the next.")
+readline()
+
+pe = heatmap(xc, yc, ue, title = "Exact Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pe)
+println("Press Enter to close the plot and exit.")
+readline()
+
+max_err = maximum(abs, ue - ua)
+println("Maximum error: $max_err")
+rel_err = 100 * max_err ./ (maximum(ue) - minimum(ue))
+println("Relative error: $rel_err")

julia/MOLE.jl/examples/elliptic2DXPerYDirichlet.jl

@@ -0,0 +1,51 @@
+using LinearAlgebra
+using SparseArrays
+using Plots
+
+import MOLE: Operators, BCs
+
+# Parameters
+k = 2        # Order of accuracy
+m = 20       # Number of cells in x-direction
+n = m + 1    # Number of cells in y-direction
+dx = 1.0 / m # Step size in x-direction
+dy = 1.0 / n # Step size in y-direction
+
+# Grid
+xc = (dx / 2.0) : dx : (1.0 - dx / 2.0)
+yc = [0; (dy / 2.0) : dy : (1.0 - dy / 2.0); 1.0]
+X = ones(1, n + 2) .* xc
+Y = yc' .* ones(m, 1)
+
+# Exact Solution
+ue = Y .* (1.0 .- Y) .* sin.(2.0 * pi .* X)
+
+# Boundary Conditions
+dc = (0.0, 0.0, 1.0, 1.0)
+nc = (0.0, 0.0, 0.0, 0.0)
+v = ([0.0], [0.0], vec(zeros(m, 1)), vec(zeros(m, 1)))
+bc = BCs.ScalarBC2D(dc, nc, v)
+
+# Operator
+A = - Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+
+# RHS
+b = 2.0 * sin.(2.0 * pi .* X) .* (1.0 .+ 2.0 * pi^2 .* Y .* (1.0 .- Y))
+b = vec(b)
+
+# Apply BCs
+A0, b0 = BCs.addScalarBC!(sparse(A), b, bc, k, m, dx, n, dy)
+
+# Approximate Solution
+ua = A0 \ b0
+ua = Matrix(reshape(ua, m, n + 2))
+
+pa = heatmap(xc, yc, ua, title = "Approximate Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pa)
+println("Press Enter to close the plot and open the next.")
+readline()
+
+pe = heatmap(xc, yc, ue, title = "Exact Solution", xlabel = "X", ylabel = "Y", colorbar_title = "u(x,y)", aspect_ratio = :equal)
+display(pe)
+println("Press Enter to close the plot and exit.")
+readline()

julia/MOLE.jl/examples/parabolic2D.jl

@@ -0,0 +1,74 @@
+using LinearAlgebra
+using SparseArrays
+using Plots
+
+import MOLE: Operators, BCs
+
+# Parameters
+k = 2               # Order of accuracy
+nu = 0.1            # Diffusion coefficient
+m = 40              # Number of cells in x-direction
+n = 50              # Number of cells in y-direction
+dx = 2.0 / m        # Step size in x-direction
+dy = 2.0 / n        # Step size in y-direction
+t = 3.0             # Final time
+method = "implicit" # Method of Euler time integration
+dt = 0.0            # Step size in time
+if method == "explicit"; dt = 0.001; else; dt = 0.01; end
+
+# Grid
+xc = [0; (dx / 2.0) : dx : (2.0 - dx / 2.0); 2.0]
+yc = [0; (dy / 2.0) : dy : (2.0 - dy / 2.0); 2.0]
+
+# Initial Conditions
+u = zeros(m + 2, n + 2)
+for i = 1:m
+    for j = 1:n
+        if ((1.0 ≤ yc[j]) && (yc[j] ≤ 1.5) && (1.0 ≤ xc[i]) && (xc[i] ≤ 1.5))
+            u[i, j] = 2.0
+        end
+    end
+end
+u = Matrix(reshape(u, :, 1))
+
+# Boundary Conditions
+dc = (1.0, 1.0, 1.0, 1.0)
+nc = (0.0, 0.0, 0.0, 0.0)
+v = (vec(zeros(n, 1)), vec(zeros(n, 1)), vec(zeros(m + 2, 1)), vec(zeros(m + 2, 1)))
+bc = BCs.ScalarBC2D(dc, nc, v)
+
+# Operator
+L = Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc)
+if method == "explicit"
+    L = nu * dt .* sparse(L) .+ sparse(Matrix(I, size(L)))
+else
+    L = sparse(Matrix(I, size(L))) .- nu * dt .* sparse(L)
+end
+
+for it = 0:(t/dt)
+    heatmap(xc, yc, u,
+        size = (800, 600),
+        xlabel = "x",
+        ylabel = "y",
+        title = "$method\n2-D Diffusion with ν = $nu -- time (t) = $(it*dt)",
+        xlims = (0, 2),
+        ylims = (0, 2),
+        clims = (0, 2),
+        colorbar = true,
+        colorbar_title = "u(x,y)"
+    )
+    display(current())
+
+    global u = reshape(u', :, 1)
+
+    if method == "explicit"
+        global u = L * u
+    else
+        global u = L \ u
+    end
+
+    global u = reshape(u, n + 2, m + 2)'
+end
+
+println("Press Enter to close the plot and exit.")
+readline()

julia/MOLE.jl/src/BCs/robinBC.jl

@@ -30,4 +30,45 @@ function robinBC(k::Int, m::Int, dx, a, b)
     G = Operators.grad(k,m,dx)
 
     BC = A + B*G;
+end
+
+"""
+    robinBC2D(k, m, dx, n, dy, a, b)
+
+Returns a two-dimensional mimetic boundary condition operator that imposes a boundary condition of Robin's type
+
+# Arguments
+- `k::Int`: Order of accuracy
+- `m::Int`: Number of cells in x-direction
+- `dx`: Step size in x-direction
+- `n::Int`: Number of cells in y-direction
+- `dy`: Step size in y-direction
+- `a`: Dirichlet Coefficient
+- `b`: Neumann Coefficient
+"""
+function robinBC2D(k::Int, m::Int, dx, n::Int, dy, a, b)
+
+    Bm = robinBC(k, m, dx, a, b)
+    Bn = robinBC(k, n, dy, a, b)
+
+    Im = Matrix(I, m + 2, m + 2)
+
+    In = Matrix(I, n + 2, n + 2)
+    In[1, 1] = 0
+    In[end, end] = 0
+
+    BC1 = kron(In, Bm)
+    BC2 = kron(Bn, Im)
+
+    BC = BC1 + BC2
+
+end
+
+"""
+    robinBC(k, m, dx, n, dy, a, b)
+
+Alias of robinBC2D
+"""
+function robinBC(k::Int, m::Int, dx, n::Int, dy, a, b)
+    return robinBC2D(k, m, dx, n, dy, a, b);
 end