Adding interpol.jl, hyperbolic1D.jl and edit on Operators.jl

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)

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)   

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
+``` 

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
+```
+

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.

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"),
+)
+

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")
+)
+
+

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
+end # module Operators

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
+