diff --git a/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/Burgers1D.md b/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/Burgers1D.md new file mode 100644 index 000000000..db8342f54 --- /dev/null +++ b/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/Burgers1D.md @@ -0,0 +1,19 @@ +# Burgers 1D + +Solves the conservative form of the inviscid Burger's Equation in 1D. + +```math +\frac{\partial }{\partial t} U + \frac{\partial}{\partial x} \left(\frac{U^{2}}{2}\right)=0 +``` + +with $U = u(x,t)$ defined on the domain $x \in [-15,15]$, from time $t \in [0,10]$ and initial condition + +```math +u(x, 0) = e^{\frac{-x^{2}}{50}} +``` + +The wave is allowed to propagate across the domain while the area under the curve is calculated. + + +--- +This example is implemented in [`burgers1D.jl`](https://github.com/csrc-sdsu/mole/blob/main/julia/MOLE/jl/examples/hyperbolic/burgers1D.jl) diff --git a/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/Hyperbolic1D.md b/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/Hyperbolic1D.md new file mode 100644 index 000000000..45603cffe --- /dev/null +++ b/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/Hyperbolic1D.md @@ -0,0 +1,36 @@ +# Hyperbolic 1D + +Solves the 1D advection equation with periodic boundary conditions. + +```math +\frac{\partial }{\partial t} U + a \frac{\partial }{\partial x} U = 0 +``` + +where $U = u(x,t)$ and $a = 1$ is the advection velocity. The domain $x \in [0,1]$ and $t \in [0,1]$ with initial condition + +```math +u(x, 0) = \sin(2\pi x) +``` + +Periodic boundary conditions are used + +```math +u(0, t) = u(1, t) +``` + +Using finite differences for the time derivative + +```math +\frac{\partial U}{\partial t} = \frac{U^{n+1}_{i}-U^{n}_{i}}{\delta t} +``` + +where $U_{i}^{n}$ is $u(x_{i},t_{n})$ and the mimetic operator $\mathbf{D}$ for the space derivative. + +```math +\frac{U^{n+1}_{i}-U^{n}_{i}}{\delta t} + a \mathbf{D}U_{i}^{n} = 0\\ +\frac{U^{n+1}_{i}-U^{n}_{i}}{\delta t} = - a \mathbf{D}U_{i}^{n}\\ +U^{n+1}_{i} = U^{n}_{i} - a \delta t \mathbf{D} U_{i}^{n}\\ +``` + +--- +This example is implemented in [`hyperbolic1D.jl`](https://github.com/csrc-sdsu/mole/blob/main/julia/MOLE.jl/examples/hyperbolic/hyperbolic1D.jl) diff --git a/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/index.md b/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/index.md new file mode 100644 index 000000000..b5bad5fb3 --- /dev/null +++ b/julia/MOLE.jl/docs/src/examples/Hyperbolic/1D/index.md @@ -0,0 +1,10 @@ +# 1D Hyperbolic Problems + +These examples demonstrate the solution of hyperbolic equations in one dimension. + +```@contents +:maxdepth: 1 + +Hyperbolic1D +Burgers1D +``` diff --git a/julia/MOLE.jl/docs/src/examples/Hyperbolic/index.md b/julia/MOLE.jl/docs/src/examples/Hyperbolic/index.md new file mode 100644 index 000000000..cfa485ed0 --- /dev/null +++ b/julia/MOLE.jl/docs/src/examples/Hyperbolic/index.md @@ -0,0 +1,11 @@ +# Hyperbolic Problems + +Hyperbolic partial differential equations model wave propagation and transport phenomena where information travels along characteristic curves. Examples include the wave equation, transport equation, and conservation laws such as Burgers' equation. + +```@content +:maxdepth: 2 +:caption: Contents + +1D/index +``` + diff --git a/julia/MOLE.jl/docs/src/examples/index.md b/julia/MOLE.jl/docs/src/examples/index.md index 7c374e54a..c89576b79 100644 --- a/julia/MOLE.jl/docs/src/examples/index.md +++ b/julia/MOLE.jl/docs/src/examples/index.md @@ -11,6 +11,7 @@ The name for OCTAVE/MATLAB, C++, and Julia will be the same (i.e., the files `el src/examples/Elliptic/index src/examples/Parabolic/index +src/examples/Hyperbolic/index ``` More examples will be added in the future. diff --git a/julia/MOLE.jl/docs/src/index.md b/julia/MOLE.jl/docs/src/index.md index 3aeef4709..c7c0d7434 100644 --- a/julia/MOLE.jl/docs/src/index.md +++ b/julia/MOLE.jl/docs/src/index.md @@ -132,6 +132,10 @@ Currently, the following examples are available in the MOLE Julia package. - 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. +- Hyperbolic Problems + - 1D Examples + - ```burgers1D```: A script that solves the 1D conservative form of inviscid Burgers equation using mimetic operators. + - ```hyperbolic1D```: A script that solves the 1D Hyperbolic Equation with Periodic boundary conditions 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. diff --git a/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl b/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl new file mode 100644 index 000000000..8e1a73551 --- /dev/null +++ b/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl @@ -0,0 +1,106 @@ +#= +Solve the 1D Inviscid Burgers' equation +Upwind Scheme is used and the equation is written in conservative form. +Initial condition: exp(-x^2/50) +=# + +using Plots +import MOLE: Operators, BCs + + +function upwind_burgers(u0, xgrid, T, k, m, dx) +#Solves Burgers equation using MOLE Operators + #= + INPUTS: + u0 : Initial condition + xgrid : staggered grid + T : final time + k : Operator order of accuarcy + m : number of cells + OUTPUTS: + U : solution + t : time interval + =# + #CFL Condition for explicit schemes + dt = dx; + + #Create stepsize for time with given T + t = collect(0.0:dt:T) + + #Use of MOLE Operators + D = Operators.div(k, m, dx) + I = Operators.interpol(m, 1.0) + + #Premultiply out of time loop (does not change) + D = -dt/2*D*I + + #Create an array that holds solution at each time step + U = zeros(length(t)+1,length(xgrid)) + + #Set initial condition at first time element + U[1,:] .= u0.(xgrid) + + for k in eachindex(t) + U[k+1,:] .= U[k,:] + D * U[k,:].^2; + end + + return t, U + +end + +### MAIN LOOP ### + +#Domain limits +west = -15; +east = 15; + +#number of cells +m = 300; + +#stepsize +dx = (east - west)/m; + +#Simulation time +T = 13 + +#Operator's order of accuracy +k = 2; + +# 1D Staggered grid +xgrid = [west; (west + (dx/2)):dx: (east - (dx/2)); east]; + +# Initial Condition +u0(x) = exp.(-(x.^2)/50); + +t, U = upwind_burgers(u0, xgrid, T, k, m, dx) + +path = joinpath(@__DIR__, "output") # Output path to store generated plots +mkpath(path) + + +animation = Plots.@animate for k in eachindex(t) + Plots.plot(xgrid, U[k,:]; + label = "approximated", + xlabel = "x", + ylabel = "u(x,t)", + grid = true, + title = "1D Inviscid Burgers' Equation t = $(round(t[k]; digits=3))", + legend = :topleft + ) +end + + +Plots.gif(animation, joinpath(path,"burgers1D.gif"), fps=10) + +Plots.png( + Plots.plot(xgrid, U[end,:]; + label = "approximated", + xlabel = "x", + ylabel = "u(x,t)", + grid = true, + title = "1D Inviscid Burgers' Equation t = $(round(t[end]; digits=3))", + legend = :topleft + ), + joinpath(path,"burgers_end.png"), +) + diff --git a/julia/MOLE.jl/examples/hyperbolic/hyperbolic1D.jl b/julia/MOLE.jl/examples/hyperbolic/hyperbolic1D.jl new file mode 100644 index 000000000..b70b87d86 --- /dev/null +++ b/julia/MOLE.jl/examples/hyperbolic/hyperbolic1D.jl @@ -0,0 +1,136 @@ +#= +Julia implementation of the 1D Advection equation with periodic boundary conditions +from MATLAB MOLE example using the Leapfrog Scheme +=# + +using Plots +import MOLE: Operators, BCs + +function advection_hyperbolic(u0, a, grid, T, k, m, dx) + #Solves the 1D Advection equation using MOLE Operators + #= + INPUTS: + u0 : Initial condition + a : velocity + grid : staggered grid + T : final time + k : Operator order of accuracy + m : number of cells + OUTPUTS: + U : solution + t : time interval + =# + #CFL condition for explicit schemes + dt = dx/abs(a); + + #Create stepsize for time with given T + t = collect(0.0:dt:T) + + #Use of MOLE Operators + D = Operators.div(k,m,dx); + I = Operators.interpol(m,0.5); + + #Periodic Boundary Conditions imposed on Divergence Operator + D[1, 2] = 1/(2*dx); + D[1, end-1] = -1/(2*dx); + D[end, 2] = 1/(2*dx); + D[end, end-1] = -1/(2*dx); + + #Premultiply out of time loop (does not change) + D = -a*dt*2 *D*I; + + #Create an array that holds solution at each time step + U = zeros(length(t)+2, length(grid)) + + #Set initial condition at first time element + U[1,:] .= u0.(grid) + + #Leapfrog scheme requires two steps, hence we would use Euler's step + U[2,:] .= U[1,:] + D/2*U[1,:]; + + for k in eachindex(t) + U[k+2,:] .= U[k,:] + D * U[k+1,:]; + end + + return t, U +end + +### MAIN LOOP ### + +#Domain limits +west = 0; +east = 1; + +#velocity +a = 1; + +#number of cells +m = 50; + +#stepsize +dx = (east - west) / m; + +#Simulation time +T = 1; + +#Operator's order of accuracy +k = 2; + +#1D Staggered grid +xgrid = [west; (west + (dx/2)):dx: (east - (dx/2)); east]; + +#Initial Condition +u0(x) = sin.(2 * π * x); + +t, U = advection_hyperbolic(u0, a, xgrid, T, k, m, dx) + +#Output path to store generated plot +path = joinpath(@__DIR__, "output") + +#Creation of animation +animation = Plots.@animate for k in eachindex(t) + Plots.plot(xgrid, U[k,:]; + label = "approximated", + xlabel = "x", + ylabel = "u(x,t)", + grid = true, + ls = :dot, + lw = 3, + title = "1D Advection Equation with Periodic BC t = $(round(t[k], sigdigits=2))", + legend = :bottomright, + xlims = (west, east), + ylims = (-1.5, 1.5) + ) + Plots.plot!(xgrid, sin.(2*π .*(xgrid .- a*t[k])); + label = "exact", + ) +end + +#Storing animation as gif to output directory +Plots.gif(animation, joinpath(path, "hyperbolic1D.gif"), fps=10) + +#Creating plot object for approximation +p1 = plot(xgrid, U[end-2,:]; #Due to length of t and U being off by 2 + label = "approximated", + xlabel = "x", + ylabel = "u(x,t)", + grid = true, + ls = :dot, + lw = 3, + title = "1D Advection Equation with Periodic BC t = $(round(t[end], sigdigits=2))", + legend = :bottomright, + xlims = (west, east), + ylims = (-1.5, 1.5) + ) +#Mutating plot to include exact solution +plot!(xgrid, sin.(2*π .*(xgrid .- a*t[end])); + label = "exact", + ) + +#Storing png file of last iteration +Plots.png( + p1, + joinpath(path, "hyperbolic_end.png") +) + + diff --git a/julia/MOLE.jl/examples/hyperbolic/output/burgers1D.gif b/julia/MOLE.jl/examples/hyperbolic/output/burgers1D.gif new file mode 100644 index 000000000..77f1515b2 Binary files /dev/null and b/julia/MOLE.jl/examples/hyperbolic/output/burgers1D.gif differ diff --git a/julia/MOLE.jl/examples/hyperbolic/output/burgers_end.png b/julia/MOLE.jl/examples/hyperbolic/output/burgers_end.png new file mode 100644 index 000000000..4bd51c61b Binary files /dev/null and b/julia/MOLE.jl/examples/hyperbolic/output/burgers_end.png differ diff --git a/julia/MOLE.jl/examples/hyperbolic/output/hyperbolic1D.gif b/julia/MOLE.jl/examples/hyperbolic/output/hyperbolic1D.gif new file mode 100644 index 000000000..ff980e4b4 Binary files /dev/null and b/julia/MOLE.jl/examples/hyperbolic/output/hyperbolic1D.gif differ diff --git a/julia/MOLE.jl/examples/hyperbolic/output/hyperbolic_end.png b/julia/MOLE.jl/examples/hyperbolic/output/hyperbolic_end.png new file mode 100644 index 000000000..a13b08ac8 Binary files /dev/null and b/julia/MOLE.jl/examples/hyperbolic/output/hyperbolic_end.png differ diff --git a/julia/MOLE.jl/src/Operators/Operators.jl b/julia/MOLE.jl/src/Operators/Operators.jl index 20491d46b..7af4cfb75 100644 --- a/julia/MOLE.jl/src/Operators/Operators.jl +++ b/julia/MOLE.jl/src/Operators/Operators.jl @@ -3,5 +3,6 @@ module Operators include("gradient.jl") include("divergence.jl") include("laplacian.jl") +include("interpol.jl") -end # module Operators \ No newline at end of file +end # module Operators diff --git a/julia/MOLE.jl/src/Operators/interpol.jl b/julia/MOLE.jl/src/Operators/interpol.jl new file mode 100644 index 000000000..9d47b74cc --- /dev/null +++ b/julia/MOLE.jl/src/Operators/interpol.jl @@ -0,0 +1,41 @@ +using SparseArrays +#= + SPDX-License-Identifier: GPL-3.0-or-later + © 2008-2024 San Diego State University Research Foundation (SDSURF). + See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details. +=# + +function interpol(m::Int,c::Float64) + """ + implementation of MATLAB interpol.m operator for Julia + INPUTS: + 'm::Int' : number of cells + 'c::Float' : left interpolation coefficient + OUTPUTS: + I : a (m+1)×(m+2) one-dimensional interpolator of 2nd-order + """ + + #Assertions: + @assert m>=4 ["m >= 4"]; + @assert c >= 0 && c <= 1 ["0 <= c <= 1"]; + + #Dimensions of I: + n_rows = m+1; + n_cols = m+2; + + I = zeros(n_rows, n_cols); + + I[1,1] = 1; + I[end,end]=1; + + #Average between two continuous cells + avg = [c 1-c]; + + j = 2; + for i in 2: n_rows - 1 + I[i, j:j+2-1] = avg; + j = j + 1; + end + return I +end + diff --git a/julia/MOLE.jl/test/Operators/interpol.jl b/julia/MOLE.jl/test/Operators/interpol.jl new file mode 100644 index 000000000..6f2742c21 --- /dev/null +++ b/julia/MOLE.jl/test/Operators/interpol.jl @@ -0,0 +1,43 @@ +tol = 1e-10 +#= +@testset "Testing interpolating operator for m=5 and c=0.5" begin + m = 5 + c = 0.5 + I = Operators.interpol(m,c) + A = [1 0 0 0 0 0 0; + 0 0.5 0.5 0 0 0 0; + 0 0 0.5 0.5 0 0 0; + 0 0 0 0.5 0.5 0 0; + 0 0 0 0 0.5 0.5 0; + 0 0 0 0 0 0 1] +@test norm(I-A) < tol +end + + +@testset "Testing interpolating operator for m=5 and c=1" begin + m = 5 + c = 1.0 + I = Operators.interpol(m,c) + A = [1 0 0 0 0 0 0; + 0 1 0 0 0 0 0; + 0 0 1 0 0 0 0; + 0 0 0 1 0 0 0; + 0 0 0 0 1 0 0; + 0 0 0 0 0 0 1] +@test norm(I-A) < tol +end +=# + +@testset "Testing interpolator operator for m=5 and c=$c" for c=0:0.25:1 + m = 5 + I = Operators.interpol(m,c) + A = [1 0 0 0 0 0 0; + 0 c 1-c 0 0 0 0; + 0 0 c 1-c 0 0 0; + 0 0 0 c 1-c 0 0; + 0 0 0 0 c 1-c 0; + 0 0 0 0 0 0 1] +@test norm(I-A) < tol +end + + diff --git a/julia/MOLE.jl/test/runtests.jl b/julia/MOLE.jl/test/runtests.jl index 037ea604e..46d705942 100644 --- a/julia/MOLE.jl/test/runtests.jl +++ b/julia/MOLE.jl/test/runtests.jl @@ -14,6 +14,10 @@ import MOLE: Operators, BCs @testset "Testing 1D laplacian" begin include("Operators/laplacian.jl") end + + @testset "Testing interpol" begin + include("Operators/interpol.jl") + end end @testset "Testing Boundary Conditions" begin