From 777db6a3d49797b0c601468a4edbbef9a0711840 Mon Sep 17 00:00:00 2001 From: Valeria Barra Date: Wed, 24 Jun 2026 15:24:40 -0700 Subject: [PATCH 1/2] MOLE.jl: Apply JuliaFormatter --- julia/MOLE.jl/docs/make.jl | 13 +- julia/MOLE.jl/examples/elliptic/elliptic1D.jl | 18 +-- .../elliptic/elliptic1DaddScalarBC.jl | 26 ++-- .../elliptic2DXDirichletYDirichlet.jl | 46 +++---- .../elliptic/elliptic2DXPerYDirichlet.jl | 44 +++---- .../MOLE.jl/examples/hyperbolic/burgers1D.jl | 73 ++++++----- .../examples/hyperbolic/hyperbolic1D.jl | 96 +++++++-------- .../MOLE.jl/examples/parabolic/parabolic2D.jl | 24 ++-- julia/MOLE.jl/src/BCs/BCs.jl | 2 +- julia/MOLE.jl/src/BCs/robinBC.jl | 16 +-- julia/MOLE.jl/src/BCs/scalarBC.jl | 108 +++++++++++----- julia/MOLE.jl/src/MOLE.jl | 2 +- julia/MOLE.jl/src/Operators/divergence.jl | 94 ++++++++------ julia/MOLE.jl/src/Operators/gradient.jl | 116 ++++++++++-------- julia/MOLE.jl/src/Operators/interpol.jl | 11 +- julia/MOLE.jl/src/Operators/laplacian.jl | 28 +++-- julia/MOLE.jl/test/BCs/scalarBC.jl | 84 +++++++------ julia/MOLE.jl/test/Operators/divergence.jl | 24 ++-- julia/MOLE.jl/test/Operators/gradient.jl | 22 ++-- julia/MOLE.jl/test/Operators/interpol.jl | 18 ++- julia/MOLE.jl/test/Operators/laplacian.jl | 18 +-- julia/MOLE.jl/test/runtests.jl | 2 +- 22 files changed, 492 insertions(+), 393 deletions(-) diff --git a/julia/MOLE.jl/docs/make.jl b/julia/MOLE.jl/docs/make.jl index 4e3d6a21c..c9873af50 100644 --- a/julia/MOLE.jl/docs/make.jl +++ b/julia/MOLE.jl/docs/make.jl @@ -1,4 +1,4 @@ -push!(LOAD_PATH,"../src/") +push!(LOAD_PATH, "../src/") using Documenter, MOLE, SparseArrays makedocs( @@ -12,15 +12,12 @@ makedocs( "1D Elliptic Problems" => [ "Overview" => "examples/Elliptic/1D/index.md", "Elliptic 1D" => "examples/Elliptic/1D/Elliptic1D.md", - "Elliptic 1D Add Scalar Boundary Conditions" => - "examples/Elliptic/1D/Elliptic1D-add-Scalar-BC.md", + "Elliptic 1D Add Scalar Boundary Conditions" => "examples/Elliptic/1D/Elliptic1D-add-Scalar-BC.md", ], "2D Elliptic Problems" => [ "Overview" => "examples/Elliptic/2D/index.md", - "X Dirichlet Y Dirichlet" => - "examples/Elliptic/2D/Elliptic2D-X-Dirichlet-Y-Dirichlet.md", - "X Periodic Y Dirichlet" => - "examples/Elliptic/2D/Elliptic2D-X-Periodic-Y-Dirichlet.md", + "X Dirichlet Y Dirichlet" => "examples/Elliptic/2D/Elliptic2D-X-Dirichlet-Y-Dirichlet.md", + "X Periodic Y Dirichlet" => "examples/Elliptic/2D/Elliptic2D-X-Periodic-Y-Dirichlet.md", ], ], "Hyperbolic Problems" => [ @@ -40,4 +37,4 @@ makedocs( ], ], ], -) \ No newline at end of file +) diff --git a/julia/MOLE.jl/examples/elliptic/elliptic1D.jl b/julia/MOLE.jl/examples/elliptic/elliptic1D.jl index 1211483e3..60629cb1a 100644 --- a/julia/MOLE.jl/examples/elliptic/elliptic1D.jl +++ b/julia/MOLE.jl/examples/elliptic/elliptic1D.jl @@ -28,15 +28,15 @@ dx = (east-west)/m # step length path = joinpath(@__DIR__, "output") # Output path to store generated plots mkpath(path) -L = Operators.lap(k,m,dx) +L = Operators.lap(k, m, dx) # Impose Robin boundary condition on laplacian operator a = 1.0 b = 1.0 -L = L + BCs.robinBC(k,m,dx,a,b) +L = L + BCs.robinBC(k, m, dx, a, b) # 1D staggered grid -grid = [west; (west+(dx/2)):dx:(east-(dx/2)); east] +grid = [west; (west + (dx / 2)):dx:(east - (dx / 2)); east] # RHS U = exp.(grid) @@ -46,13 +46,13 @@ U[end] = 2*exp(1) U = L\U # Plot result -p = Plots.scatter(grid, U, label="Approximated", show = false) -plot!(p, grid, exp.(grid), label="Analytical", show = false) +p = Plots.scatter(grid, U, label = "Approximated", show = false) +plot!(p, grid, exp.(grid), label = "Analytical", show = false) plot!( p, - xlabel="x", - ylabel="u(x)", - title="Poisson's equation with Robin BC", - show = false + xlabel = "x", + ylabel = "u(x)", + title = "Poisson's equation with Robin BC", + show = false, ) Plots.png(p, joinpath(path, "elliptic1D_output.png")) diff --git a/julia/MOLE.jl/examples/elliptic/elliptic1DaddScalarBC.jl b/julia/MOLE.jl/examples/elliptic/elliptic1DaddScalarBC.jl index 026061db9..19437ee27 100644 --- a/julia/MOLE.jl/examples/elliptic/elliptic1DaddScalarBC.jl +++ b/julia/MOLE.jl/examples/elliptic/elliptic1DaddScalarBC.jl @@ -24,15 +24,15 @@ import MOLE: Operators, BCs west = 0.0 east = 1.0 -k = 6 # operator order of accuracy -m = 2k + 1 # minimum number of cells to attain desired accuracy +k = 6 # operator order of accuracy +m = 2k + 1 # minimum number of cells to attain desired accuracy dx = (east - west) / m # step length path = joinpath(@__DIR__, "output") # Output path to store generated plots mkpath(path) # 1D staggered grid -grid = [west; (west + dx/2):dx:(east - dx/2); east] +grid = [west; (west + dx / 2):dx:(east - dx / 2); east] # Mimetic Laplacian operator L = Operators.lap(k, m, dx) @@ -40,7 +40,7 @@ L = Operators.lap(k, m, dx) # Robin boundary conditions: dc*u + nc*(du/dn) = v dc = (1.0, 1.0) nc = (1.0, 1.0) -v = (0.0, 2 * exp(1.0)) +v = (0.0, 2 * exp(1.0)) bc = BCs.ScalarBC1D(dc, nc, v) @@ -54,14 +54,14 @@ L0, rhs0 = BCs.addScalarBC!(sparse(L), rhs, bc, k, m, dx) u = L0 \ rhs0 # Plot -p = scatter(grid, u; label="Approximated", marker=:circle, show = false) -plot!(p, grid, exp.(grid); label="Analytical", show = false) +p = scatter(grid, u; label = "Approximated", marker = :circle, show = false) +plot!(p, grid, exp.(grid); label = "Analytical", show = false) plot!( - p, - title="Poisson's equation with Robin BC using addScalarBC", - xlabel="x", - ylabel="u(x)", - legend=:topleft, - show = false + p, + title = "Poisson's equation with Robin BC using addScalarBC", + xlabel = "x", + ylabel = "u(x)", + legend = :topleft, + show = false, ) -Plots.png(p, joinpath(path, "elliptic1DaddScalarBC_output.png")) \ No newline at end of file +Plots.png(p, joinpath(path, "elliptic1DaddScalarBC_output.png")) diff --git a/julia/MOLE.jl/examples/elliptic/elliptic2DXDirichletYDirichlet.jl b/julia/MOLE.jl/examples/elliptic/elliptic2DXDirichletYDirichlet.jl index 5d8de2782..dcef9b1b1 100644 --- a/julia/MOLE.jl/examples/elliptic/elliptic2DXDirichletYDirichlet.jl +++ b/julia/MOLE.jl/examples/elliptic/elliptic2DXDirichletYDirichlet.jl @@ -34,8 +34,8 @@ path = joinpath(@__DIR__, "output") # Output path to store generated plots mkpath(path) # Grid -xc = [0; (dx / 2.0) : dx : (pi - dx / 2.0); pi] -yc = [0; (dy / 2.0) : dy : (pi - dy / 2.0); pi] +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) @@ -45,11 +45,11 @@ 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])) +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) +A = - Operators.lap(k, m, dx, n, dy, dc = dc, nc = nc) # RHS b = zeros(m + 2, n + 2) @@ -64,37 +64,37 @@ ua = Matrix(reshape(ua, m + 2, n + 2)) Plots.png( Plots.heatmap( - xc, - yc, - ua, - title = "Approximate Solution", - xlabel = "X", - ylabel = "Y", - colorbar_title = "u(x,y)", + xc, + yc, + ua, + title = "Approximate Solution", + xlabel = "X", + ylabel = "Y", + colorbar_title = "u(x,y)", aspect_ratio = :equal, colormap = :jet1, - show = false + show = false, ), - joinpath(path, "elliptic2DXDYD_approximate.png") + joinpath(path, "elliptic2DXDYD_approximate.png"), ) Plots.png( Plots.heatmap( - xc, - yc, - ue, - title = "Exact Solution", - xlabel = "X", - ylabel = "Y", - colorbar_title = "u(x,y)", + xc, + yc, + ue, + title = "Exact Solution", + xlabel = "X", + ylabel = "Y", + colorbar_title = "u(x,y)", aspect_ratio = :equal, colormap = :jet1, - show = false + show = false, ), - joinpath(path, "elliptic2DXDYD_exact.png") + joinpath(path, "elliptic2DXDYD_exact.png"), ) 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") \ No newline at end of file +println("Relative error: $rel_err") diff --git a/julia/MOLE.jl/examples/elliptic/elliptic2DXPerYDirichlet.jl b/julia/MOLE.jl/examples/elliptic/elliptic2DXPerYDirichlet.jl index 7b25f2955..66ce09704 100644 --- a/julia/MOLE.jl/examples/elliptic/elliptic2DXPerYDirichlet.jl +++ b/julia/MOLE.jl/examples/elliptic/elliptic2DXPerYDirichlet.jl @@ -34,8 +34,8 @@ path = joinpath(@__DIR__, "output") # Output path to store generated plots mkpath(path) # Grid -xc = (dx / 2.0) : dx : (1.0 - dx / 2.0) -yc = [0; (dy / 2.0) : dy : (1.0 - dy / 2.0); 1.0] +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) @@ -49,7 +49,7 @@ 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) +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)) @@ -64,37 +64,37 @@ ua = Matrix(reshape(ua, m, n + 2)) Plots.png( Plots.heatmap( - xc, - yc, - ua, - title = "Approximate Solution", - xlabel = "X", - ylabel = "Y", - colorbar_title = "u(x,y)", + xc, + yc, + ua, + title = "Approximate Solution", + xlabel = "X", + ylabel = "Y", + colorbar_title = "u(x,y)", aspect_ratio = :equal, colormap = :jet1, - show = false + show = false, ), - joinpath(path, "elliptic2DXPYD_approximate.png") + joinpath(path, "elliptic2DXPYD_approximate.png"), ) Plots.png( Plots.heatmap( - xc, - yc, - ue, - title = "Exact Solution", - xlabel = "X", - ylabel = "Y", - colorbar_title = "u(x,y)", + xc, + yc, + ue, + title = "Exact Solution", + xlabel = "X", + ylabel = "Y", + colorbar_title = "u(x,y)", aspect_ratio = :equal, colormap = :jet1, - show = false + show = false, ), - joinpath(path, "elliptic2DXPYD_exact.png") + joinpath(path, "elliptic2DXPYD_exact.png"), ) max_err = maximum(abs, ue - ua) println("Maximum error: $max_err") rel_err = 100 * max_err ./ (maximum(ue) - minimum(ue)) -println("Realative error: $rel_err") \ No newline at end of file +println("Realative error: $rel_err") diff --git a/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl b/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl index 8e1a73551..dc6d3f3ac 100644 --- a/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl +++ b/julia/MOLE.jl/examples/hyperbolic/burgers1D.jl @@ -9,7 +9,7 @@ import MOLE: Operators, BCs function upwind_burgers(u0, xgrid, T, k, m, dx) -#Solves Burgers equation using MOLE Operators + #Solves Burgers equation using MOLE Operators #= INPUTS: u0 : Initial condition @@ -22,26 +22,26 @@ function upwind_burgers(u0, xgrid, T, k, m, dx) t : time interval =# #CFL Condition for explicit schemes - dt = dx; - + dt = dx; + #Create stepsize for time with given T - t = collect(0.0:dt:T) - + t = collect(0.0:dt:T) + #Use of MOLE Operators - D = Operators.div(k, m, dx) - I = Operators.interpol(m, 1.0) - + 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 - + D = -dt/2*D*I + #Create an array that holds solution at each time step - U = zeros(length(t)+1,length(xgrid)) - + U = zeros(length(t)+1, length(xgrid)) + #Set initial condition at first time element - U[1,:] .= u0.(xgrid) - + U[1, :] .= u0.(xgrid) + for k in eachindex(t) - U[k+1,:] .= U[k,:] + D * U[k,:].^2; + U[k + 1, :] .= U[k, :] + D * U[k, :] .^ 2; end return t, U @@ -67,11 +67,11 @@ T = 13 k = 2; # 1D Staggered grid -xgrid = [west; (west + (dx/2)):dx: (east - (dx/2)); east]; +xgrid = [west; (west + (dx / 2)):dx:(east - (dx / 2)); east]; # Initial Condition -u0(x) = exp.(-(x.^2)/50); - +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 @@ -79,28 +79,27 @@ 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 - ) + 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.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"), + 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 index b70b87d86..b10bc1497 100644 --- a/julia/MOLE.jl/examples/hyperbolic/hyperbolic1D.jl +++ b/julia/MOLE.jl/examples/hyperbolic/hyperbolic1D.jl @@ -24,32 +24,32 @@ function advection_hyperbolic(u0, a, grid, T, k, m, dx) dt = dx/abs(a); #Create stepsize for time with given T - t = collect(0.0:dt:T) + t = collect(0.0:dt:T) #Use of MOLE Operators - D = Operators.div(k,m,dx); - I = Operators.interpol(m,0.5); - + 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); + 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; - + 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) + U[1, :] .= u0.(grid) #Leapfrog scheme requires two steps, hence we would use Euler's step - U[2,:] .= U[1,:] + D/2*U[1,:]; + U[2, :] .= U[1, :] + D/2*U[1, :]; for k in eachindex(t) - U[k+2,:] .= U[k,:] + D * U[k+1,:]; + U[k + 2, :] .= U[k, :] + D * U[k + 1, :]; end return t, U @@ -77,7 +77,7 @@ T = 1; k = 2; #1D Staggered grid -xgrid = [west; (west + (dx/2)):dx: (east - (dx/2)); east]; +xgrid = [west; (west + (dx / 2)):dx:(east - (dx / 2)); east]; #Initial Condition u0(x) = sin.(2 * π * x); @@ -85,52 +85,50 @@ 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") +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", - ) + 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) +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) - ) +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", - ) +plot!(xgrid, sin.(2*π .* (xgrid .- a*t[end])); + label = "exact", +) #Storing png file of last iteration Plots.png( p1, - joinpath(path, "hyperbolic_end.png") + joinpath(path, "hyperbolic_end.png"), ) - - diff --git a/julia/MOLE.jl/examples/parabolic/parabolic2D.jl b/julia/MOLE.jl/examples/parabolic/parabolic2D.jl index 1cd8069c1..9a5a03b60 100644 --- a/julia/MOLE.jl/examples/parabolic/parabolic2D.jl +++ b/julia/MOLE.jl/examples/parabolic/parabolic2D.jl @@ -48,13 +48,13 @@ path = joinpath(@__DIR__, "output") # Output path to store generated plots mkpath(path) # 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] +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 +for i in 1:m + for j in 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 @@ -66,14 +66,14 @@ u = vec(reshape(u, :, 1)) 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)) - ) + 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) +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 @@ -86,7 +86,7 @@ function time_step(L, u, t, dt, method) num_it = ceil(Int, t / dt) sol = zeros(length(u), 1 + num_it) - for it = 0:num_it + for it in 0:num_it sol[:, it + 1] = u if method == "explicit" @@ -114,7 +114,7 @@ anim = Plots.@animate for i in 1:length(sol[1, :]) colorbar = true, colorbar_title = "u(x,y)", colormap = :jet1, - show = false + show = false, ) end -Plots.gif(anim, joinpath(path, "parabolic2D.gif"), fps = ceil(Int, 1 / dt)) \ No newline at end of file +Plots.gif(anim, joinpath(path, "parabolic2D.gif"), fps = ceil(Int, 1 / dt)) diff --git a/julia/MOLE.jl/src/BCs/BCs.jl b/julia/MOLE.jl/src/BCs/BCs.jl index a27d8ca55..c35d80d98 100644 --- a/julia/MOLE.jl/src/BCs/BCs.jl +++ b/julia/MOLE.jl/src/BCs/BCs.jl @@ -5,4 +5,4 @@ using SparseArrays include("robinBC.jl") include("scalarBC.jl") -end # module BCs \ No newline at end of file +end # module BCs diff --git a/julia/MOLE.jl/src/BCs/robinBC.jl b/julia/MOLE.jl/src/BCs/robinBC.jl index 9bd266cd1..7ac6a0f41 100644 --- a/julia/MOLE.jl/src/BCs/robinBC.jl +++ b/julia/MOLE.jl/src/BCs/robinBC.jl @@ -19,15 +19,15 @@ Returns a m+2 by m+2 one-dimensional mimetic boundary operator that imposes a bo - `b`: Neumann Coefficient """ function robinBC(k::Int, m::Int, dx, a, b) - A = zeros(m+2,m+2) - A[1,1] = a - A[m+2,m+2] = a + A = zeros(m+2, m+2) + A[1, 1] = a + A[m + 2, m + 2] = a - B = zeros(m+2,m+1) - B[1,1] = -b - B[m+2,m+1] = b + B = zeros(m+2, m+1) + B[1, 1] = -b + B[m + 2, m + 1] = b - G = Operators.grad(k,m,dx) + G = Operators.grad(k, m, dx) BC = A + B*G; end @@ -71,4 +71,4 @@ 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 \ No newline at end of file +end diff --git a/julia/MOLE.jl/src/BCs/scalarBC.jl b/julia/MOLE.jl/src/BCs/scalarBC.jl index 77f44f6dc..8542cbf63 100644 --- a/julia/MOLE.jl/src/BCs/scalarBC.jl +++ b/julia/MOLE.jl/src/BCs/scalarBC.jl @@ -24,9 +24,9 @@ abstract type AbstractScalarBC{D} end - v: prescribed boundary value g (left,right) """ struct ScalarBC1D{T} <: AbstractScalarBC{1} - dc::NTuple{2,T} - nc::NTuple{2,T} - v::NTuple{2,T} + dc::NTuple{2, T} + nc::NTuple{2, T} + v::NTuple{2, T} end @@ -39,9 +39,9 @@ end - vc: prescribed boundary value g (left, right, bottom, top) """ struct ScalarBC2D{T} <: AbstractScalarBC{2} - dc::NTuple{4,T} - nc::NTuple{4,T} - v::NTuple{4,AbstractVector{T}} + dc::NTuple{4, T} + nc::NTuple{4, T} + v::NTuple{4, AbstractVector{T}} end @@ -53,14 +53,26 @@ end Internal helper: build LHS contributions for 1D BC. Returns (Al, Ar). """ -function _scalarbc1d_lhs(k::Integer, m::Integer, dx, dc::NTuple{2,T}, nc::NTuple{2,T}) where {T} +function _scalarbc1d_lhs( + k::Integer, + m::Integer, + dx, + dc::NTuple{2, T}, + nc::NTuple{2, T}, +) where {T} n = m + 2 Al = spzeros(T, n, n) Ar = spzeros(T, n, n) - if dc[1] != zero(T); Al[1, 1] = dc[1]; end - if dc[2] != zero(T); Ar[end, end] = dc[2]; end + if dc[1] != zero(T) + ; + Al[1, 1] = dc[1]; + end + if dc[2] != zero(T) + ; + Ar[end, end] = dc[2]; + end Bl = spzeros(T, n, m + 1) Br = spzeros(T, n, m + 1) @@ -68,8 +80,14 @@ function _scalarbc1d_lhs(k::Integer, m::Integer, dx, dc::NTuple{2,T}, nc::NTuple Gl = grad(k, m, dx) Gr = grad(k, m, dx) - if nc[1] != zero(T); Bl[1, 1] = -nc[1]; end - if nc[2] != zero(T); Br[end, end] = nc[2]; end + if nc[1] != zero(T) + ; + Bl[1, 1] = -nc[1]; + end + if nc[2] != zero(T) + ; + Br[end, end] = nc[2]; + end return Al + Bl * Gl, Ar + Br * Gr end @@ -77,7 +95,7 @@ end """ Internal helper: overwrite RHS at boundary indices. """ -@inline function _scalarbc_rhs!(b, vec::NTuple{2,Int}, v::NTuple{2,T}) where {T} +@inline function _scalarbc_rhs!(b, vec::NTuple{2, Int}, v::NTuple{2, T}) where {T} b[vec[1]] = v[1] b[vec[2]] = v[2] return b @@ -88,12 +106,13 @@ end Signature keeps the discretization params (`k,m,dx`) separate from `bc`. """ function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC1D{T}, - k::Integer, m::Integer, dx) where {T} + k::Integer, m::Integer, dx) where {T} dc, nc, v = bc.dc, bc.nc, bc.v # Equivalent of MATLAB: q = find(dc.^2 + nc.^2, 1) - hasbc = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) + hasbc = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) if !hasbc return A, b end @@ -104,7 +123,13 @@ function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC1D{T}, # We remove the existing nonzeros in those rows. sub = A[[vec[1], vec[2]], :] rows, cols, vals = findnz(sub) # rows are 1..2 in the submatrix - A .-= sparse((rows .== 1) .* vec[1] .+ (rows .== 2) .* vec[2], cols, vals, size(A,1), size(A,2)) + A .-= sparse( + (rows .== 1) .* vec[1] .+ (rows .== 2) .* vec[2], + cols, + vals, + size(A, 1), + size(A, 2), + ) # Zero out boundary entries of b b[vec[1]] = zero(eltype(b)) @@ -129,12 +154,22 @@ end Internal helper: build LHS contributions for 2D BC. Returns (Al, Ar, Ab, At). """ -function _scalarbc2D_lhs(k::Integer, m::Integer, dx, n::Integer, dy, dc::NTuple{4,T}, nc::NTuple{4,T}) where {T} +function _scalarbc2D_lhs( + k::Integer, + m::Integer, + dx, + n::Integer, + dy, + dc::NTuple{4, T}, + nc::NTuple{4, T}, +) where {T} Abcl, Abcr, Abcb, Abct = 0, 0, 0, 0 # periodic case - hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) - hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) + hasbclr = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) + hasbcbt = + (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) if hasbclr @@ -177,10 +212,21 @@ end """ Internal helper: overwrite RHS at boundary indices """ -@inline function _scalarbc2d_rhs(b, dc::NTuple{4,T}, nc::NTuple{4,T}, v::NTuple{4,AbstractVector{T}}, rl::AbstractVector{Int}, rr::AbstractVector{Int}, rb::AbstractVector{Int}, rt::AbstractVector{Int}) where {T} - - hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) - hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) +@inline function _scalarbc2d_rhs( + b, + dc::NTuple{4, T}, + nc::NTuple{4, T}, + v::NTuple{4, AbstractVector{T}}, + rl::AbstractVector{Int}, + rr::AbstractVector{Int}, + rb::AbstractVector{Int}, + rt::AbstractVector{Int}, +) where {T} + + hasbclr = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) + hasbcbt = + (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) if hasbclr b[rl] = v[1] @@ -200,12 +246,14 @@ end Signature keeps the discretization params (`k,m,dx,n,dy`) separate from `bc`. """ function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC2D{T}, - k::Integer, m::Integer, dx, n::Integer, dy) where {T} + k::Integer, m::Integer, dx, n::Integer, dy) where {T} dc, nc, v = bc.dc, bc.nc, bc.v - hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) - hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) + hasbclr = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) + hasbcbt = + (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) rl, rr, rb, rt = [0], [0], [0], [0] # periodic case @@ -222,8 +270,8 @@ function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC2D{T}, Abc1 = Abcl .+ Abcr rowsbc1, _, _ = findnz(Abc1) # For some reason, the below code doesn't work: - # rows1, cols1, s1 = findnz(A[rowsbc1, :]) - # A = A .- sparse(rows1, cols1, s1, size(A, 1), size(A, 2)) + # rows1, cols1, s1 = findnz(A[rowsbc1, :]) + # A = A .- sparse(rows1, cols1, s1, size(A, 1), size(A, 2)) A[rowsbc1, :] = spzeros(size(A[rowsbc1, :])) # This does though # Update matrix A with boundary information A = A .+ Abc1 @@ -233,7 +281,7 @@ function addScalarBC!(A::SparseMatrixCSC, b::AbstractVector, bc::ScalarBC2D{T}, end if hasbcbt - + rb, _, _ = findnz(Abcb) rt, _, _ = findnz(Abct) rb = unique(rb) @@ -264,7 +312,7 @@ end Alias for `addScalarBC!(A, b, bc, k, m, dx, n, dy)` """ function addScalarBC2D!(A::SparseMatrixCSC, b::AbstractVector{T}, bc::ScalarBC2D{T}, - k::Integer, m::Integer, dx, n::Integer, dy) where {T} + k::Integer, m::Integer, dx, n::Integer, dy) where {T} return addScalarBC!(A, b, bc, k, m, dx, n, dy) end @@ -277,4 +325,4 @@ end Convenience constructor from vectors/tuples. Accepts 2-element vectors too. """ -ScalarBC1D(dc, nc, v) = ScalarBC1D(tuple(dc...), tuple(nc...), tuple(v...)) \ No newline at end of file +ScalarBC1D(dc, nc, v) = ScalarBC1D(tuple(dc...), tuple(nc...), tuple(v...)) diff --git a/julia/MOLE.jl/src/MOLE.jl b/julia/MOLE.jl/src/MOLE.jl index 5360c41cf..77ad9d077 100644 --- a/julia/MOLE.jl/src/MOLE.jl +++ b/julia/MOLE.jl/src/MOLE.jl @@ -9,4 +9,4 @@ module MOLE include("Operators/Operators.jl") include("BCs/BCs.jl") -end # module MOLE \ No newline at end of file +end # module MOLE diff --git a/julia/MOLE.jl/src/Operators/divergence.jl b/julia/MOLE.jl/src/Operators/divergence.jl index 7d5af46f7..1ce1dd1d5 100644 --- a/julia/MOLE.jl/src/Operators/divergence.jl +++ b/julia/MOLE.jl/src/Operators/divergence.jl @@ -20,9 +20,16 @@ Returns a one-dimensional mimetic divergence operator. Default is non periodic. - `dc::NTuple{2,T}`: Dirichlet coefficients of left and right boundaries (optional) - `nc::NTuple{2,T}`: Neumann coefficients of left and right boundaries (optional) """ -function div(k::Int, m::Int, dx; dc::NTuple{2,T} = (1.0, 1.0), nc::NTuple{2,T} = (1.0, 1.0)) where {T} - - hasbc = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) +function div( + k::Int, + m::Int, + dx; + dc::NTuple{2, T} = (1.0, 1.0), + nc::NTuple{2, T} = (1.0, 1.0), +) where {T} + + hasbc = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) if hasbc D = divNonPeriodic(k, m, dx) return D @@ -69,35 +76,38 @@ function divNonPeriodic(k::Int, m::Int, dx) throw(DomainError(m, "m must be >= 2*k + 1")) end - D = zeros(m+2,m+1) - if k == 2 - for i = 2:(m+1) - D[i,(i-1):i] = [-1 1] + D = zeros(m+2, m+1) + if k == 2 + for i in 2:(m + 1) + D[i, (i - 1):i] = [-1 1] end elseif k == 4 A = [-11/12 17/24 3/8 -5/24 1/24] - D[2,1:5] = A - D[m+1, (m-3):(m+1)] = -reverse(A) - for i = 3:m - D[i,(i-2):(i+1)] = [1/24 -9/8 9/8 -1/24] + D[2, 1:5] = A + D[m + 1, (m - 3):(m + 1)] = -reverse(A) + for i in 3:m + D[i, (i - 2):(i + 1)] = [1/24 -9/8 9/8 -1/24] end elseif k == 6 - A = [-1627/1920 211/640 59/48 -235/192 91/128 -443/1920 31/960; - 31/960 -687/640 129/128 19/192 -3/32 21/640 -3/640] - D[2:3,1:7] = A - D[m:(m+1), (m-5):(m+1)] = -rot180(A) - for i = 4:(m-1) - D[i,(i-3):(i+2)] = [-3/640 25/384 -75/64 75/64 -25/384 3/640] + A = [-1627/1920 211/640 59/48 -235/192 91/128 -443/1920 31/960; + 31/960 -687/640 129/128 19/192 -3/32 21/640 -3/640] + D[2:3, 1:7] = A + D[m:(m + 1), (m - 5):(m + 1)] = -rot180(A) + for i in 4:(m - 1) + D[i, (i - 3):(i + 2)] = [-3/640 25/384 -75/64 75/64 -25/384 3/640] end elseif k == 8 - A = [-1423/1792 -491/7168 7753/3072 -18509/5120 3535/1024 -2279/1024 953/1024 -1637/7168 2689/107520; - 2689/107520 -36527/35840 4259/5120 6497/15360 -475/1024 1541/5120 -639/5120 1087/35840 -59/17920; - -59/17920 1175/21504 -1165/1024 1135/1024 25/3072 -251/5120 25/1024 -45/7168 5/7168] - D[2:4, 1:9] = A; - D[2:4,1:9] = A - D[(m-1):(m+1), (m-7):(m+1)] = -rot180(A) - for i = 5:(m-2) - D[i,(i-4):(i+3)] = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168] + A = [ + -1423/1792 -491/7168 7753/3072 -18509/5120 3535/1024 -2279/1024 953/1024 -1637/7168 2689/107520; + 2689/107520 -36527/35840 4259/5120 6497/15360 -475/1024 1541/5120 -639/5120 1087/35840 -59/17920; + -59/17920 1175/21504 -1165/1024 1135/1024 25/3072 -251/5120 25/1024 -45/7168 5/7168 + ] + D[2:4, 1:9] = A; + D[2:4, 1:9] = A + D[(m - 1):(m + 1), (m - 7):(m + 1)] = -rot180(A) + for i in 5:(m - 2) + D[i, (i - 4):(i + 3)] = + [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168] end end D = (1/dx)*D; @@ -117,7 +127,7 @@ Returns a m by m periodic mimetic divergence operator. function divPeriodic(k::Int, m::Int, dx) D = - gradPeriodic(k, m, dx)'; - + end @@ -137,9 +147,9 @@ function divNonUniform(k::Int, ticks::AbstractVector) m, _ = size(D) if size(ticks, 1) == 1 - J = diagm((D * ticks').^-1) + J = diagm((D * ticks') .^ -1) else - J = diagm((D * ticks).^-1) + J = diagm((D * ticks) .^ -1) end D = J * D @@ -167,15 +177,25 @@ Returns a two-dimensional mimetic divergence operator. Default is non periodic. - `dc::NTuple{4,T}`: Dirichlet coefficients of left, right, bottom, and top boundaries (optional) - `nc::NTuple{4,T}`: Neumann coefficients of left, right, bottom, and top boundaries (optional) """ -function div(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0), nc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0)) where {T} - - hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) - hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) +function div( + k::Int, + m::Int, + dx, + n::Int, + dy; + dc::NTuple{4, T} = (1.0, 1.0, 1.0, 1.0), + nc::NTuple{4, T} = (1.0, 1.0, 1.0, 1.0), +) where {T} + + hasbclr = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) + hasbcbt = + (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) if hasbclr Dx = divNonPeriodic(k, m, dx) Im = Matrix(I, m + 2, m + 2) - Im = Im[:, 2:end-1] + Im = Im[:, 2:(end - 1)] else Dx = divPeriodic(k, m, dx) Im = Matrix(I, m, m) @@ -184,7 +204,7 @@ function div(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T} = (1.0, 1.0, 1.0, 1 if hasbcbt Dy = divNonPeriodic(k, n, dy) In = Matrix(I, n + 2, n + 2) - In = In[:, 2:end-1] + In = In[:, 2:(end - 1)] else Dy = divPeriodic(k, n, dy) In = Matrix(I, n, n) @@ -233,12 +253,12 @@ function div2DNonUniform(k::Int, xticks::AbstractVector, yticks::AbstractVector) Im = Matrix(I, m, m) In = Matrix(I, n, n) - Im = Im[:, 2:end-1] - In = In[:, 2:end-1] + Im = Im[:, 2:(end - 1)] + In = In[:, 2:(end - 1)] Sx = kron(In, Dx) Sy = kron(Dy, Im) D = [Sx Sy] -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/src/Operators/gradient.jl b/julia/MOLE.jl/src/Operators/gradient.jl index 0f6a2a211..05286ccd5 100644 --- a/julia/MOLE.jl/src/Operators/gradient.jl +++ b/julia/MOLE.jl/src/Operators/gradient.jl @@ -22,9 +22,16 @@ Returns a one-dimensional mimetic gradient operator. Default is non periodic. - `dc::NTuple{2,T}`: Dirichlet coefficients of the left and right boundaries (optional) - `nc::NTuple{2,T}`: Neumann coefficients of the left and right boundaries (optional) """ -function grad(k::Int, m::Int, dx; dc::NTuple{2,T} = (1.0, 1.0), nc::NTuple{2,T} = (1.0, 1.0)) where {T} - - hasbc = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) +function grad( + k::Int, + m::Int, + dx; + dc::NTuple{2, T} = (1.0, 1.0), + nc::NTuple{2, T} = (1.0, 1.0), +) where {T} + + hasbc = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) if hasbc G = gradNonPeriodic(k, m, dx) return G @@ -73,40 +80,43 @@ function gradNonPeriodic(k::Int, m::Int, dx) throw(DomainError(m, "m must be >= 2*k")) end - G = zeros(m+1,m+2) + G = zeros(m+1, m+2) if k == 2 A = [-8/3 3 -1/3] - G[1,1:3] = A #[-8/3 3 -1/3] - G[m+1,m:(m+2)] = -reverse(A) #[1/3, -3, 8/3] - for i = 2:m - G[i,i:(i+1)] = [-1 1] + G[1, 1:3] = A #[-8/3 3 -1/3] + G[m + 1, m:(m + 2)] = -reverse(A) #[1/3, -3, 8/3] + for i in 2:m + G[i, i:(i + 1)] = [-1 1] end elseif k == 4 - A = [-352/105 35/8 -35/24 21/40 -5/56; - 16/105 -31/24 29/24 -3/40 1/168] - G[1:2,1:5] = A - G[m:(m+1), (m-2):(m+2)] = -rot180(A) - for i = 3:(m-1) - G[i,(i-1):(i+2)] = [1/24 -9/8 9/8 -1/24] + A = [-352/105 35/8 -35/24 21/40 -5/56; + 16/105 -31/24 29/24 -3/40 1/168] + G[1:2, 1:5] = A + G[m:(m + 1), (m - 2):(m + 2)] = -rot180(A) + for i in 3:(m - 1) + G[i, (i - 1):(i + 2)] = [1/24 -9/8 9/8 -1/24] end elseif k == 6 - A = [-13016/3465 693/128 -385/128 693/320 -495/448 385/1152 -63/1408; - 496/3465 -811/640 449/384 -29/960 -11/448 13/1152 -37/21120; - -8/385 179/1920 -153/128 381/320 -101/1344 1/128 -3/7040]; - G[1:3,1:7] = A - G[(m-1):(m+1), (m-4):(m+2)] = -rot180(A) - for i = 4:(m-2) - G[i,(i-2):(i+3)] = [-3/640 25/384 -75/64 75/64 -25/384 3/640] + A = [-13016/3465 693/128 -385/128 693/320 -495/448 385/1152 -63/1408; + 496/3465 -811/640 449/384 -29/960 -11/448 13/1152 -37/21120; + -8/385 179/1920 -153/128 381/320 -101/1344 1/128 -3/7040]; + G[1:3, 1:7] = A + G[(m - 1):(m + 1), (m - 4):(m + 2)] = -rot180(A) + for i in 4:(m - 2) + G[i, (i - 2):(i + 3)] = [-3/640 25/384 -75/64 75/64 -25/384 3/640] end elseif k == 8 - A = [-182144/45045 6435/1024 -5005/1024 27027/5120 -32175/7168 25025/9216 -12285/11264 3465/13312 -143/5120; - 86048/675675 -131093/107520 49087/46080 10973/76800 -4597/21504 4019/27648 -10331/168960 2983/199680 -2621/1612800; - -3776/225225 8707/107520 -17947/15360 29319/25600 -533/21504 -263/9216 903/56320 -283/66560 257/537600; - 32/9009 -543/35840 265/3072 -1233/1024 8625/7168 -775/9216 639/56320 -15/13312 1/21504] - G[1:4,1:9] = A - G[(m-2):(m+1), (m-6):(m+2)] = -rot180(A) - for i = 5:(m-3) - G[i,(i-3):(i+4)] = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168] + A = [ + -182144/45045 6435/1024 -5005/1024 27027/5120 -32175/7168 25025/9216 -12285/11264 3465/13312 -143/5120; + 86048/675675 -131093/107520 49087/46080 10973/76800 -4597/21504 4019/27648 -10331/168960 2983/199680 -2621/1612800; + -3776/225225 8707/107520 -17947/15360 29319/25600 -533/21504 -263/9216 903/56320 -283/66560 257/537600; + 32/9009 -543/35840 265/3072 -1233/1024 8625/7168 -775/9216 639/56320 -15/13312 1/21504 + ] + G[1:4, 1:9] = A + G[(m - 2):(m + 1), (m - 6):(m + 2)] = -rot180(A) + for i in 5:(m - 3) + G[i, (i - 3):(i + 4)] = + [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168] end end G = (1/dx)*G; @@ -124,8 +134,8 @@ Returns a m by m periodic mimetic gradient operator - `dx`: Step size in x-direction """ function gradPeriodic(k::Int, m::Int, dx) - - if k < 2 || k > 8 + + if k < 2 || k > 8 throw(DomainError(k, "k must be >= 2 and <= 8")) end @@ -136,11 +146,11 @@ function gradPeriodic(k::Int, m::Int, dx) if m < 2*k throw(DomainError(m, "m must be >= 2*k")) end - + V = zeros(1, m) idx = fill(-1, m, m) idx[:, 1] = 1:m - idx = cumsum(idx, dims=2) + idx = cumsum(idx, dims = 2) idx = mod.(idx .+ m, m) .+ 1 if k == 2 @@ -159,7 +169,7 @@ function gradPeriodic(k::Int, m::Int, dx) elseif k == 8 V[1, 1:6] = [-245/3072, 1225/1024, -1225/1024, 245/3072, -49/5120, 5/7168] - V[1, end-1:end] = [-5/7168, 49/5120] + V[1, (end - 1):end] = [-5/7168, 49/5120] end @@ -183,9 +193,9 @@ function gradNonUniform(k::Int, ticks::AbstractVector) G = grad(k, length(ticks) - 2, 1) if size(ticks, 1) == 1 - J = diagm((G*ticks').^-1) + J = diagm((G*ticks') .^ -1) else - J = diagm((G*ticks).^-1) + J = diagm((G*ticks) .^ -1) end G = J * G @@ -210,15 +220,25 @@ Returns a two-dimensional mimetic gradient operator. Default is non periodic. - `dc::NTuple{4,T}`: Dirichlet coefficients of the left, right, bottom, and top boundaries (optional) - `nc::NTuple{4,T}`: Neumann coefficients of the left, right, bottom, and top boundaries (optional) """ -function grad(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0), nc::NTuple{4,T} = (1.0, 1.0, 1.0, 1.0)) where {T} - - hasbclr = (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) - hasbcbt = (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) +function grad( + k::Int, + m::Int, + dx, + n::Int, + dy; + dc::NTuple{4, T} = (1.0, 1.0, 1.0, 1.0), + nc::NTuple{4, T} = (1.0, 1.0, 1.0, 1.0), +) where {T} + + hasbclr = + (dc[1] != zero(T)) || (dc[2] != zero(T)) || (nc[1] != zero(T)) || (nc[2] != zero(T)) + hasbcbt = + (dc[3] != zero(T)) || (dc[4] != zero(T)) || (nc[3] != zero(T)) || (nc[4] != zero(T)) if hasbclr Gx = gradNonPeriodic(k, m, dx) Im = Matrix(I, m + 2, m + 2) - Im = Im[2:end-1, :] + Im = Im[2:(end - 1), :] else Gx = gradPeriodic(k, m, dx) Im = Matrix(I, m, m) @@ -227,7 +247,7 @@ function grad(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4,T} = (1.0, 1.0, 1.0, if hasbcbt Gy = gradNonPeriodic(k, n, dy) In = Matrix(I, n + 2, n + 2) - In = In[2:end-1, :] + In = In[2:(end - 1), :] else Gy = gradPeriodic(k, n, dy) In = Matrix(I, n, n) @@ -276,14 +296,14 @@ function gradNonPeriodic(k::Int, m::Int, dx, n::Int, dy) Im = Matrix(I, m + 2, m + 2) In = Matrix(I, n + 2, n + 2) - Im = Im[2:end-1, :] - In = In[2:end-1, :] + Im = Im[2:(end - 1), :] + In = In[2:(end - 1), :] Sx = kron(In, Gx) Sy = kron(Gy, Im) G = [Sx; Sy]; - + end @@ -336,12 +356,12 @@ function gradNonUniform(k::Int, xticks::AbstractVector, yticks::AbstractVector) Im = Matrix(I, m, m) In = Matrix(I, n, n) - Im = Im[2:end-1, :] - In = In[2:end-1, :] + Im = Im[2:(end - 1), :] + In = In[2:(end - 1), :] Sx = kron(In, Gx) Sy = kron(Gy, Im) G = [Sx; Sy] -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/src/Operators/interpol.jl b/julia/MOLE.jl/src/Operators/interpol.jl index 99c0ad804..1b96a5c9d 100644 --- a/julia/MOLE.jl/src/Operators/interpol.jl +++ b/julia/MOLE.jl/src/Operators/interpol.jl @@ -16,7 +16,7 @@ Returns a (m+1)×(m+2) one-dimensional interpolator of 2nd-order - `c::Float`: left interpolation coefficient """ function interpol(m::Int, c::Float64) - + #Assertions: @assert m>=4 ["m >= 4"]; @assert c >= 0 && c <= 1 ["0 <= c <= 1"]; @@ -27,17 +27,16 @@ function interpol(m::Int, c::Float64) I = zeros(n_rows, n_cols); - I[1,1] = 1; - I[end,end]=1; + 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; + 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/src/Operators/laplacian.jl b/julia/MOLE.jl/src/Operators/laplacian.jl index 6aed4c083..1cd8c71a8 100644 --- a/julia/MOLE.jl/src/Operators/laplacian.jl +++ b/julia/MOLE.jl/src/Operators/laplacian.jl @@ -20,9 +20,15 @@ Returns a m+2 by m+2 one-dimensional mimetic laplacian operator. Default is non - `dc::NTuple{2,T}`: Dirichlet coefficients of the left and right boundaries (optional) - `nc::NTuple{2,T}`: Neumann coefficients of the left and right boundaries (optional) """ -function lap(k::Int, m::Int, dx; dc::NTuple{2,T} = (1.0, 1.0), nc::NTuple{2,T} = (1.0, 1.0)) where {T} - D = div(k, m, dx, dc=dc, nc=nc) - G = grad(k, m, dx, dc=dc, nc=nc) +function lap( + k::Int, + m::Int, + dx; + dc::NTuple{2, T} = (1.0, 1.0), + nc::NTuple{2, T} = (1.0, 1.0), +) where {T} + D = div(k, m, dx, dc = dc, nc = nc) + G = grad(k, m, dx, dc = dc, nc = nc) L = D*G; end @@ -46,11 +52,19 @@ Returns a two-dimensional mimetic laplacian operator. Default is non periodic. - `dc::NTuple{4,T}`: Dirichlet coefficients of the left and right boundaries (optional) - `nc::NTuple{4,T}`: Neumann coefficients of the left and right boundaries (optional) """ -function lap(k::Int, m::Int, dx, n::Int, dy; dc::NTuple{4}=(1.0, 1.0, 1.0, 1.0), nc::NTuple{4}=(1.0, 1.0, 1.0, 1.0)) +function lap( + k::Int, + m::Int, + dx, + n::Int, + dy; + dc::NTuple{4} = (1.0, 1.0, 1.0, 1.0), + nc::NTuple{4} = (1.0, 1.0, 1.0, 1.0), +) - D = div(k, m, dx, n, dy, dc=dc, nc=nc) - G = grad(k, m, dx, n, dy, dc=dc, nc=nc) + D = div(k, m, dx, n, dy, dc = dc, nc = nc) + G = grad(k, m, dx, n, dy, dc = dc, nc = nc) L = D*G -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/test/BCs/scalarBC.jl b/julia/MOLE.jl/test/BCs/scalarBC.jl index f167c5be3..dc604e369 100644 --- a/julia/MOLE.jl/test/BCs/scalarBC.jl +++ b/julia/MOLE.jl/test/BCs/scalarBC.jl @@ -2,19 +2,19 @@ @testset "no BCs test" begin # Problem size - m = 8 - n = m + 2 - k = 2 + m = 8 + n = m + 2 + k = 2 dx = 0.1 # Nontrivial sparse A and b - A_dense = reshape(collect(1.0:(n*n)), n, n) ./ 13.0 + A_dense = reshape(collect(1.0:(n * n)), n, n) ./ 13.0 A = sparse(A_dense) b = collect(1.0:n) ./ 7.0 dc0 = (0.0, 0.0) nc0 = (0.0, 0.0) - v0 = (9.0, 9.0) # required by ScalarBC1D even if unused + v0 = (9.0, 9.0) # required by ScalarBC1D even if unused bc0 = BCs.ScalarBC1D(dc0, nc0, v0) A2, b2 = BCs.addScalarBC!(A, b, bc0, k, m, dx) @@ -26,19 +26,19 @@ @testset "Dirichlet/Neumann BC with grad and dense matrix reference" begin # Problem size - m = 8 - n = m + 2 - k = 2 + m = 8 + n = m + 2 + k = 2 dx = 0.1 # Nontrivial sparse A and b - A_dense = reshape(collect(1.0:(n*n)), n, n) ./ 13.0 + A_dense = reshape(collect(1.0:(n * n)), n, n) ./ 13.0 A = sparse(A_dense) b = collect(1.0:n) ./ 7.0 dc = (2.0, 3.0) nc = (4.0, 5.0) - v = (7.0, 8.0) + v = (7.0, 8.0) bc = BCs.ScalarBC1D(dc, nc, v) A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx) @@ -64,7 +64,7 @@ Bl = zeros(Float64, n, m + 1) Br = zeros(Float64, n, m + 1) Bl[1, 1] = -nc[1] - Br[end, end] = nc[2] + Br[end, end] = nc[2] A_ref .+= (Al + Bl * Gl) .+ (Ar + Br * Gr) @@ -72,25 +72,25 @@ b_ref[1] = v[1] b_ref[end] = v[2] - @test isapprox(norm(A2 - sparse(A_ref)), 0.0; rtol=1e-12, atol=1e-12) - @test isapprox(norm(b2 - b_ref), 0.0; rtol=1e-12, atol=1e-12) + @test isapprox(norm(A2 - sparse(A_ref)), 0.0; rtol = 1e-12, atol = 1e-12) + @test isapprox(norm(b2 - b_ref), 0.0; rtol = 1e-12, atol = 1e-12) end @testset "Dirichlet/Neumann BC with sparse format reference" begin # Problem size - m = 8 - n = m + 2 - k = 2 + m = 8 + n = m + 2 + k = 2 dx = 0.1 # Nontrivial sparse A and b - A_dense = reshape(collect(1.0:(n*n)), n, n) ./ 13.0 + A_dense = reshape(collect(1.0:(n * n)), n, n) ./ 13.0 A = sparse(A_dense) b = collect(1.0:n) ./ 7.0 dc = (2.0, 3.0) nc = (4.0, 5.0) - v = (7.0, 8.0) + v = (7.0, 8.0) bc = BCs.ScalarBC1D(dc, nc, v) A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx) @@ -130,7 +130,7 @@ Bl = spzeros(Float64, n, m + 1) Br = spzeros(Float64, n, m + 1) Bl[1, 1] = -nc[1] - Br[end, end] = nc[2] + Br[end, end] = nc[2] A_ref .+= (Al + Bl * Gl) .+ (Ar + Br * Gr) @@ -161,7 +161,7 @@ end # Non trivial sparse A and b A = sparse(rand(m * n, m * n)) - b = rand(m * n,) + b = rand(m * n) A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx, n, dy) @@ -183,13 +183,13 @@ end # Boundary conditions dc = (3.0, 2.0, 1.0, 4.0) nc = (1.0, 8.0, 3.0, 6.0) - v = (rand(n,), rand(n,), rand(m + 2,), rand(m + 2,)) + v = (rand(n), rand(n), rand(m + 2), rand(m + 2)) bc = BCs.ScalarBC2D(dc, nc, v) # Non trivial sparse A and b A_dense = rand(numC, numC) A = sparse(A_dense) - b = rand(numC,) + b = rand(numC) A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx, n, dy) @@ -229,7 +229,7 @@ end Inn = copy(In) Inn[1, 1] = 0 Inn[end, end] = 0 - + Bl = zeros(Float64, m + 2, m + 1) Br = zeros(Float64, m + 2, m + 1) Bb = zeros(Float64, n + 2, n + 1) @@ -252,17 +252,20 @@ end A_ref .+= Al .+ Ar .+ Ab .+ At - bdry_v = zeros(size(bdry_rows,)) + bdry_v = zeros(size(bdry_rows)) bdry_v[1:length(v[3])] = v[3] midd = [] - for i = 1:length(v[1]); midd = [midd; v[1][i]; v[2][i]]; end - bdry_v[length(v[3])+1:length([v[1];v[2];v[3]])] = midd - bdry_v[end-length(v[4])+1:end] = v[4] + for i in 1:length(v[1]) + ; + midd = [midd; v[1][i]; v[2][i]]; + end + bdry_v[(length(v[3]) + 1):length([v[1]; v[2]; v[3]])] = midd + bdry_v[(end - length(v[4]) + 1):end] = v[4] b_ref[bdry_rows] = bdry_v - @test isapprox(norm(A2 - sparse(A_ref)), 0.0; rtol=1e-12, atol=1e-12) - @test isapprox(norm(b2 - b_ref), 0.0; rtol=1e-12, atol=1e-12) + @test isapprox(norm(A2 - sparse(A_ref)), 0.0; rtol = 1e-12, atol = 1e-12) + @test isapprox(norm(b2 - b_ref), 0.0; rtol = 1e-12, atol = 1e-12) end @@ -278,12 +281,12 @@ end # Nontrivial sparse A and b A_dense = rand(numC, numC) A = sparse(A_dense) - b = rand(numC,) + b = rand(numC) # Boundary conditions dc = (1.0, 3.3, 4.0, 2.0) nc = (1.0, 3.3, 4.0, 2.0) - v = (rand(n,), rand(n,), rand(m + 2,), rand(m + 2,)) + v = (rand(n), rand(n), rand(m + 2), rand(m + 2)) bc = BCs.ScalarBC2D(dc, nc, v) A2, b2 = BCs.addScalarBC!(A, b, bc, k, m, dx, n, dy) @@ -303,17 +306,20 @@ end bdry_y[1, 1] = 1 bdry_y[end, end] = 1 bdry_rows = findall(vec(any(kron(In, bdry_x) .+ kron(bdry_y, Im) .!= 0, dims = 2))) - + A_ref[bdry_rows, :] = spzeros(size(A_ref[bdry_rows, :])) # RHS boundaries b_ref = copy(b) - bdry_v = zeros(size(bdry_rows,)) + bdry_v = zeros(size(bdry_rows)) bdry_v[1:length(v[3])] = v[3] midd = [] - for i = 1:length(v[1]); midd = [midd; v[1][i]; v[2][i]]; end - bdry_v[length(v[3])+1:length([v[1];v[2];v[3]])] = midd - bdry_v[end-length(v[4])+1:end] = v[4] + for i in 1:length(v[1]) + ; + midd = [midd; v[1][i]; v[2][i]]; + end + bdry_v[(length(v[3]) + 1):length([v[1]; v[2]; v[3]])] = midd + bdry_v[(end - length(v[4]) + 1):end] = v[4] b_ref[bdry_rows] = bdry_v # LHS boundaries @@ -335,7 +341,7 @@ end Inn = copy(In) Inn[1, 1] = 0 Inn[end, end] = 0 - + Bl = spzeros(Float64, m + 2, m + 1) Br = spzeros(Float64, m + 2, m + 1) Bb = spzeros(Float64, n + 2, n + 1) @@ -356,11 +362,11 @@ end Ab = kron(Ab, Im) At = kron(At, Im) - A_ref .+= Al .+ Ar .+ Ab .+ At + A_ref .+= Al .+ Ar .+ Ab .+ At # pure sparse comparisons @test norm(A2 - A_ref) ≤ 1e-12 @test norm(b2 - b_ref) ≤ 1e-12 end -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/test/Operators/divergence.jl b/julia/MOLE.jl/test/Operators/divergence.jl index 2bc7aee86..abcfc9046 100644 --- a/julia/MOLE.jl/test/Operators/divergence.jl +++ b/julia/MOLE.jl/test/Operators/divergence.jl @@ -1,25 +1,25 @@ tol = 1e-10 -@testset "Testing non periodic 1-D divergence for order k=$k" for k=2:2:8 +@testset "Testing non periodic 1-D divergence for order k=$k" for k in 2:2:8 m = 2*k+1 D = Operators.div(k, m, 1/m) - field = ones(m+1,1) + field = ones(m+1, 1) sol = D*field @test norm(sol) < tol end -@testset "Testing periodic 1-D divergence for order k=$k" for k=2:2:8 +@testset "Testing periodic 1-D divergence for order k=$k" for k in 2:2:8 m = 2 * k + 1 dx = 1.0 / (m - 1) dc = (0.0, 0.0) nc = (0.0, 0.0) - D = Operators.div(k, m, dx, dc=dc, nc=nc) + D = Operators.div(k, m, dx, dc = dc, nc = nc) field = ones(m, 1) sol = D * field @test norm(sol) < tol end -@testset "Testing non uniform 1-D divergence for order k=$k" for k=2:2:8 +@testset "Testing non uniform 1-D divergence for order k=$k" for k in 2:2:8 m = 2 * k + 1 ticks = sort(rand(m + 1)) D = Operators.div(k, ticks) @@ -28,7 +28,7 @@ end @test norm(sol) < tol end -@testset "Testing non periodic 2-D divergence for order k=$k" for k=2:2:8 +@testset "Testing non periodic 2-D divergence for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m @@ -39,33 +39,33 @@ end @test norm(sol) < tol end -@testset "Testing periodic 2-D divergence for order k=$k" for k=2:2:8 +@testset "Testing periodic 2-D divergence for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / (m - 1) dy = 1.0 / (n - 1) dc = (0.0, 0.0, 0.0, 0.0) nc = (0.0, 0.0, 0.0, 0.0) - D = Operators.div(k, m, dx, n, dy, dc=dc, nc=nc) + D = Operators.div(k, m, dx, n, dy, dc = dc, nc = nc) field = ones(2*m*n, 1) sol = D * field @test norm(sol) < tol end -@testset "Testing mixed periodicity 2-D divergence for order k=$k" for k=2:2:8 +@testset "Testing mixed periodicity 2-D divergence for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / (m - 1) dy = 1.0 / n dc = (0.0, 0.0, 1.0, 1.0) nc = (0.0, 0.0, 0.0, 0.0) - D = Operators.div(k, m, dx, n, dy, dc=dc, nc=nc) + D = Operators.div(k, m, dx, n, dy, dc = dc, nc = nc) field = ones(m*n + m*(n+1), 1) sol = D * field @test norm(sol) < tol end -@testset "Testing non uniform 2-D divergence for order k=$k" for k=2:2:8 +@testset "Testing non uniform 2-D divergence for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 xticks = sort(rand(m + 1)) @@ -74,4 +74,4 @@ end field = ones((m + 1) * n + m * (n + 1), 1) sol = D * field @test norm(sol) < tol -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/test/Operators/gradient.jl b/julia/MOLE.jl/test/Operators/gradient.jl index 5f52d92af..2e180b70a 100644 --- a/julia/MOLE.jl/test/Operators/gradient.jl +++ b/julia/MOLE.jl/test/Operators/gradient.jl @@ -1,6 +1,6 @@ tol = 1e-10 -@testset "Testing non periodic 1-D gradient for order k=$k" for k=2:2:8 +@testset "Testing non periodic 1-D gradient for order k=$k" for k in 2:2:8 m = 2*k+1 G = Operators.grad(k, m, 1/m) field = ones(m+2, 1) @@ -8,18 +8,18 @@ tol = 1e-10 @test norm(sol) < tol end -@testset "Testing periodic 1-D gradient for order k=$k" for k=2:2:8 +@testset "Testing periodic 1-D gradient for order k=$k" for k in 2:2:8 m = 2 * k + 1 dx = 1.0 / (m - 1) dc = (0.0, 0.0) nc = (0.0, 0.0) - G = Operators.grad(k, m, dx, dc=dc, nc=dc) + G = Operators.grad(k, m, dx, dc = dc, nc = dc) field = ones(m, 1) sol = G * field @test norm(sol) < tol end -@testset "Testing non uniform 1-D gradient for order k=$k" for k=2:2:8 +@testset "Testing non uniform 1-D gradient for order k=$k" for k in 2:2:8 m = 2 * k + 1 ticks = sort(rand(m + 2)) G = Operators.grad(k, ticks) @@ -28,7 +28,7 @@ end @test norm(sol) < tol end -@testset "Testing non periodic 2-D gradient for order k=$k" for k=2:2:8 +@testset "Testing non periodic 2-D gradient for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m @@ -39,33 +39,33 @@ end @test norm(sol) < tol end -@testset "Testing periodic 2-D gradient for order k=$k" for k=2:2:8 +@testset "Testing periodic 2-D gradient for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m dy = 1.0 / n dc = (0.0, 0.0, 0.0, 0.0) nc = (0.0, 0.0, 0.0, 0.0) - G = Operators.grad(k, m, dx, n, dy, dc=dc, nc=nc) + G = Operators.grad(k, m, dx, n, dy, dc = dc, nc = nc) field = ones(m * n, 1) sol = G * field @test norm(sol) < tol end -@testset "Testing mixed periodicity 2-D gradient for order k=$k" for k=2:2:8 +@testset "Testing mixed periodicity 2-D gradient for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m dy = 1.0 / n dc = (0.0, 0.0, 1.0, 1.0) nc = (0.0, 0.0, 1.0, 1.0) - G = Operators.grad(k, m, dx, n, dy, dc=dc, nc=nc) + G = Operators.grad(k, m, dx, n, dy, dc = dc, nc = nc) field = ones(m * (n + 2), 1) sol = G * field @test norm(sol) < tol end -@testset "Testing non uniform 2-D gradient for order k=$k" for k=2:2:8 +@testset "Testing non uniform 2-D gradient for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 xticks = sort(rand(m + 2)) @@ -74,4 +74,4 @@ end field = ones((m + 2) * (n + 2), 1) sol = G * field @test norm(sol) < tol -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/test/Operators/interpol.jl b/julia/MOLE.jl/test/Operators/interpol.jl index 79cf0719a..e8486e8a9 100644 --- a/julia/MOLE.jl/test/Operators/interpol.jl +++ b/julia/MOLE.jl/test/Operators/interpol.jl @@ -1,15 +1,13 @@ tol = 1e-10 -@testset "Testing interpolator operator for m=5 and c=$c" for c=0:0.25:1 +@testset "Testing interpolator operator for m=5 and c=$c" for c in 0:0.25:1 m = 5 - I = Operators.interpol(m,c) + 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 + 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/Operators/laplacian.jl b/julia/MOLE.jl/test/Operators/laplacian.jl index b1d2e77fa..4b4db351e 100644 --- a/julia/MOLE.jl/test/Operators/laplacian.jl +++ b/julia/MOLE.jl/test/Operators/laplacian.jl @@ -1,6 +1,6 @@ tol = 1e-10 -@testset "Testing non periodic 1-D laplacian for order k=$k" for k=2:2:8 +@testset "Testing non periodic 1-D laplacian for order k=$k" for k in 2:2:8 m = 2*k+1 L = Operators.lap(k, m, 1/m) field = ones(m+2, 1) @@ -8,18 +8,18 @@ tol = 1e-10 @test norm(sol) < tol end -@testset "Testing periodic 1-D laplacian for order k=$k" for k=2:2:8 +@testset "Testing periodic 1-D laplacian for order k=$k" for k in 2:2:8 m = 2 * k + 1 dx = 1.0 / m dc = (0.0, 0.0) nc = (0.0, 0.0) - L = Operators.lap(k, m, dx, dc=dc, nc=nc) + L = Operators.lap(k, m, dx, dc = dc, nc = nc) field = ones(m, 1) sol = L * field @test norm(sol) < tol end -@testset "Testing non periodic 2-D laplacian for order k=$k" for k=2:2:8 +@testset "Testing non periodic 2-D laplacian for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m @@ -30,28 +30,28 @@ end @test norm(sol) < tol end -@testset "Testing periodic 2-D laplacian for order k=$k" for k=2:2:8 +@testset "Testing periodic 2-D laplacian for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m dy = 1.0 / n dc = (0.0, 0.0, 0.0, 0.0) nc = (0.0, 0.0, 0.0, 0.0) - L = Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc) + L = Operators.lap(k, m, dx, n, dy, dc = dc, nc = nc) field = ones(m * n, 1) sol = L * field @test norm(sol) < tol end -@testset "Testing mixed periodicity 2-D laplacian for order k=$k" for k=2:2:8 +@testset "Testing mixed periodicity 2-D laplacian for order k=$k" for k in 2:2:8 m = 2 * k + 1 n = m + 1 dx = 1.0 / m dy = 1.0 / n dc = (0.0, 0.0, 1.0, 1.0) nc = (0.0, 0.0, 0.0, 0.0) - L = Operators.lap(k, m, dx, n, dy, dc=dc, nc=nc) + L = Operators.lap(k, m, dx, n, dy, dc = dc, nc = nc) field = ones(m * (n + 2), 1) sol = L * field @test norm(sol) < tol -end \ No newline at end of file +end diff --git a/julia/MOLE.jl/test/runtests.jl b/julia/MOLE.jl/test/runtests.jl index 46d705942..9dbd3da58 100644 --- a/julia/MOLE.jl/test/runtests.jl +++ b/julia/MOLE.jl/test/runtests.jl @@ -21,7 +21,7 @@ import MOLE: Operators, BCs end @testset "Testing Boundary Conditions" begin - + @testset "Test addScalarBC" begin include("BCs/scalarBC.jl") end From 5566ec67690caddd30dbe7d88b0487ec33db8044 Mon Sep 17 00:00:00 2001 From: Valeria Barra Date: Wed, 24 Jun 2026 15:30:14 -0700 Subject: [PATCH 2/2] julia/MOLE.jl/README.md: Update instructions on how to run the JuliaFormatter locally --- julia/MOLE.jl/README.md | 38 ++++++++++++++++++++++++++++++++++---- 1 file changed, 34 insertions(+), 4 deletions(-) diff --git a/julia/MOLE.jl/README.md b/julia/MOLE.jl/README.md index 513478789..14a2ec108 100644 --- a/julia/MOLE.jl/README.md +++ b/julia/MOLE.jl/README.md @@ -85,8 +85,38 @@ Once you have built the documentation (either from the REPL or the command line) The MOLE library contains examples demonstrating how to use the operators, in a broad range of partial differential equations (PDEs). More information on the mathematical content can be found in the [main MOLE documentation](https://mole-docs.readthedocs.io/en/main/examples/index.html). -Currently, the following examples are available in the MOLE Julia package. +The examples are organized by PDE type and dimension in the base MOLE documentation: -- Elliptic Problems - - 1D Examples - - `elliptic1D`: A script that solves the 1D Poisson's equation with Robin boundary conditions using mimetic operators. \ No newline at end of file +- **Base Examples Documentation:** https://csrc-sdsu.github.io/mole/stable/examples/ + +Currently, not all of them are available in the MOLE Julia package. They will be added in the near future. + +## Formatter + +To test the JuliaFormatter locally, please make sure to have a working julia installation and add the `JuliaFormatter` package to your either global or local environment. At a global, system level, this would be: + +``` +julia -e 'using Pkg; Pkg.add("JuliaFormatter")' +``` + +Now that your environment has the JuliaFormatter package, do the following to apply format style changes to the MOLE.jl directory: + +1. From the `mole/julia/MOLE.jl` directory, run + +``` +julia -e 'using JuliaFormatter; format(".", overwrite=true, verbose=true) || error("Formatting issues detected")' +``` + +2. Review the generated style changes applied by the formatter with + +``` +git status +``` + +and + +``` +git diff +``` + +Stage the modified files and commit the necessary changes. \ No newline at end of file