add basic files

Project.toml

@@ -3,6 +3,15 @@ uuid = "ae0c53cd-8ee1-4e9a-acfb-47eefb3fcea7"
 authors = ["Stamatis Vretinaris"]
 version = "1.0.0-DEV"
 
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
+YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+
 [compat]
 julia = "1.10"
 

src/BoundaryConditions.jl

@@ -0,0 +1,29 @@
+module BoundaryConditions
+
+using LinearAlgebra
+using StaticArrays
+
+## This matrix is setup to implement absorbing/radiative BCs
+@inline function RR(n::AbstractVector)
+    @fastmath @inbounds begin
+        mat = @SMatrix [1/2 -n[1]/2 -n[2]/2 -n[3]/2
+                        -n[1]/2 n[1]^2/2+n[2]^2+n[3]^2 -n[1] * n[2]/2 -n[1] * n[3]/2
+                        -n[2]/2 -n[1] * n[2]/2 n[1]^2+n[3]^2+n[2]^2/2 -n[2] * n[3]/2
+                        -n[3]/2 -n[1] * n[3]/2 -n[2] * n[3]/2 n[1]^2+n[2]^2+n[3]^2/2]
+    end
+    return mat
+end
+
+@inline function apply_boundary_conditions(dU, i, j, k, N)
+    n = normal_vector(i, j, k, N)
+    return RR(n) * dU
+end
+
+function normal_vector(i, j, k, N)
+    vec = @SVector [i == 1 ? -1 : i == N ? 1 : 0,
+                    j == 1 ? -1 : j == N ? 1 : 0,
+                    k == 1 ? -1 : k == N ? 1 : 0]
+    return vec / norm(vec)
+end
+
+end

src/Grids.jl

@@ -0,0 +1,78 @@
+module Grids
+
+using StaticArrays
+import Base.length
+
+export UniformGrid, spacing, coords
+
+abstract type AbstractGrid{T} end
+
+struct UniformGrid1d{T<:Real} <: AbstractGrid{T}
+    domain::SVector{2,T}
+    ncells::Int
+end
+
+function UniformGrid1d(d::AbstractVector{T}, ncells::Int) where {T<:Real}
+    return UniformGrid1d(SVector{2,T}(d), ncells)
+end
+
+function spacing(xi::T, xn::T, ncells::Integer)::Real where {T<:Real}
+    return (xn - xi) / ncells  # (n+1 -1)
+end
+
+function spacing(x::AbstractVector{<:Real})
+    spacing(x[begin], x[end], length(x) - 1)
+end
+
+function spacing(g::UniformGrid1d{T}) where {T<:Real}
+    spacing(g.domain[1], g.domain[2], g.ncells)
+end
+
+function coords(ui::Real, uf::Real, ncells::Real)
+    du = spacing(ui, uf, ncells)
+    ns = 0:1:ncells
+    u = @. ui + ns * du
+    return collect(u)
+end
+
+function coords(domain::AbstractVector{T}, ncells::Int) where {T<:Real}
+    length(domain) === 2 ||
+        throw(DimensionMismatch("domain needs to be Vector of length 2"))
+    return coords(domain[1], domain[2], ncells)
+end
+
+function coords(g::UniformGrid1d{<:Real})
+    return coords(g.domain, g.ncells)
+end
+
+struct UniformGrid{T<:Real} <: AbstractGrid{T}
+    domain::AbstractVector{<:AbstractVector{T}}
+    ncells::AbstractVector{<:Int}
+end
+
+function coords(g::UniformGrid{<:Real})
+    return coords.(g.domain, g.ncells)
+end
+
+function UniformGrid(domains::AbstractVector{<:AbstractVector{T}},
+                     ncells::Int) where {T<:Real}
+    n = length(domains)
+    ncells2 = ncells * @SVector ones(Int64, n)
+    domains = SVector{n,SVector{2,T}}(domains)
+    return UniformGrid(domains, ncells2)
+end
+
+function UniformGrid(domain::AbstractVector{T}, ncells::Int, dim::Int) where {T<:Real}
+    ncells2 = ncells * @SVector ones(Int64, dim)
+    domains = Array{SVector{2,T}}(undef, dim)
+    for i in 1:dim
+        domains[i] = domain
+    end
+    return UniformGrid(SVector{dim,SVector{2,T}}(domains), ncells2)
+end
+
+function UniformGrid(domain::AbstractVector{T}, ncells::Int) where {T<:Real}
+    return UniformGrid1d(domain, ncells)
+end
+
+end # end of module

src/InitialData.jl

@@ -0,0 +1,50 @@
+module InitialData
+
+# 3D wave
+@inline function Φ(t::Real, x::Real, y::Real, z::Real, params)::Real
+    A = params[1]
+    n = params[3]
+    L = params[4]
+    @fastmath nl = n / L
+    return @fastmath A * sinpi(nl * x) * sinpi(nl * y) * sinpi(nl * z) * cospi(√3 * nl * t)
+end
+
+@inline function Π(t::Real, x::Real, y::Real, z::Real, params)::Real
+    A = params[1]
+    n = params[3]
+    L = params[4]
+    @fastmath nl = n / L
+    @fastmath ω = √3 * nl * π
+    return @fastmath -A * ω * sinpi(nl * x) * sinpi(nl * y) * sinpi(nl * z) *
+                     sinpi(√3 * nl * t)
+end
+
+@inline function Ψx(t::Real, x::Real, y::Real, z::Real, params)::Real
+    A = params[1]
+    n = params[3]
+    L = params[4]
+    @fastmath nl = n / L
+    return @fastmath A * nl * π * cospi(nl * x) * sinpi(nl * y) * sinpi(nl * z) *
+                     cospi(√3 * nl * t)
+end
+
+@inline function Ψy(t::Real, x::Real, y::Real, z::Real, params)::Real
+    A = params[1]
+    n = params[3]
+    L = params[4]
+    @fastmath nl = n / L
+    return @fastmath A * nl * π * cospi(nl * y) * sinpi(nl * x) * sinpi(nl * z) *
+                     cospi(√3 * nl * t)
+end
+
+@inline function Ψz(t::Real, x::Real, y::Real, z::Real, params)::Real
+    A = params[1]
+    n = params[3]
+    L = params[4]
+    @fastmath nl = n / L
+    return @fastmath A * nl * π * cospi(nl * z) * sinpi(nl * y) * sinpi(nl * x) *
+                     cospi(√3 * nl * t)
+end
+
+
+end

src/InputOutput.jl

@@ -0,0 +1,86 @@
+module InputOutput
+
+using WriteVTK
+using ReadVTK
+using YAML
+using ..Grids
+
+function get_time_grid_from_yaml(metadata_file)
+    file = YAML.load_file(metadata_file)
+    dom = file["Time"]["domain"]
+    ncells = file["Time"]["ncells"]
+    every = file["Output"]["save every"]
+    time = UniformGrid(dom, div(ncells, every))
+    return time
+end
+
+function writevtk(statevector, simpath, xcoord, ti, i, pvd)
+    vtkpath = string(simpath, "/timestep_", i)
+
+    vtk = vtk_grid(vtkpath, xcoord, xcoord, xcoord)
+
+    vtk["Phi"] = statevector[:, :, :, 1]
+    vtk["Pi"] = statevector[:, :, :, 2]
+    vtk["Psix"] = statevector[:, :, :, 3]
+    vtk["Psiy"] = statevector[:, :, :, 4]
+    vtk["Psiz"] = statevector[:, :, :, 5]
+    vtk["time"] = ti
+
+    collection_add_timestep(pvd, vtk, ti)
+    vtk_save(vtk)
+    return nothing
+end
+
+function writevtk_initialdata(statevector, simpath, xcoord)
+    pvdpath = string(simpath, "/full_simulation")
+    vtkpath = string(simpath, "/timestep_0")
+
+    pvd = paraview_collection(pvdpath)
+    vtk = vtk_grid(vtkpath, xcoord, xcoord, xcoord)
+
+    vtk["Phi"] = statevector[:, :, :, 1]
+    vtk["Pi"] = statevector[:, :, :, 2]
+    vtk["Psix"] = statevector[:, :, :, 3]
+    vtk["Psiy"] = statevector[:, :, :, 4]
+    vtk["Psiz"] = statevector[:, :, :, 5]
+    vtk["time"] = 0.0
+
+    collection_add_timestep(pvd, vtk, 0.0)
+    vtk_save(vtk)
+    return pvd
+end
+
+function save_pvd(pvd)
+    vtk_save(pvd)
+    return nothing
+end
+
+# format YAML
+function write_metadata(md_path, params, save_every, simdir)
+    metadata_filename = string(md_path, "/metadata.yaml")
+    data = "Grid:
+    domain:       [$(params.grid.domain[1]), $(params.grid.domain[2])]
+    ncells:       $(params.grid.ncells)
+    grid:         $(params.grid.ncells + 1)×$(params.grid.ncells + 1)×$(params.grid.ncells + 1)
+    total points: $((params.grid.ncells + 1)^3)
+    dx=dy=dz:     $(params.h)
+
+Time:
+    CFL:          $(params.dt/params.h)
+    domain:       [0.0, $(params.ti[end])]
+    ncells:       $(length(params.ti)-1)
+    dt:           $(params.dt)
+
+Output:
+    job ID:        $(params.jobid)
+    output path:   $simdir
+    save every:    $save_every
+\n
+"
+    println(data)
+    open(metadata_filename, "w") do file
+        write(file, data)
+    end
+end
+
+end # end of module

src/Integrator.jl

@@ -0,0 +1,105 @@
+module Integrator
+
+using ..InputOutput
+using LoopVectorization
+
+const proj_path = pkgdir(Integrator)
+if occursin("cn", gethostname()) || occursin("mn", gethostname())
+    const output_path = "/gpfs/svretinaris/ScalarWave/runs/"
+else
+    const output_path = string(proj_path, "/output/")
+end
+
+function base_path(folder)
+    return string(output_path, folder)
+end
+
+function sims_path(folder)
+    return string(output_path, folder, "/sims/")
+end
+
+@inline function RK4(rhs!::F, dt::Real, reg1, reg2, reg3, reg4, params, t::Real,
+                     indices::CartesianIndices) where {F}
+    rhs!(reg4, reg1, params, t)
+    @inbounds @fastmath @avx for i in indices
+        reg3[i] = dt * reg4[i]
+        reg1[i] = reg1[i] + reg3[i] / 2
+        reg2[i] = reg3[i]
+    end
+    rhs!(reg4, reg1, params, t)
+    @inbounds @fastmath @avx for i in indices
+        reg3[i] = dt * reg4[i]
+        reg1[i] = reg1[i] + (reg3[i] - reg2[i]) / 2
+    end
+    rhs!(reg4, reg1, params, t)
+    @inbounds @fastmath @avx for i in indices
+        reg3[i] = dt * reg4[i] - reg3[i] / 2
+        reg1[i] = reg1[i] + reg3[i]
+        reg2[i] = reg2[i] / 6 - reg3[i]
+    end
+    rhs!(reg4, reg1, params, t)
+    @inbounds @fastmath @avx for i in indices
+        reg3[i] = dt * reg4[i] + reg3[i] + reg3[i]
+        reg1[i] = reg1[i] + reg2[i] + reg3[i] / 6
+    end
+    return nothing
+end
+
+function solve(rhs, statevector, params, t, grid, save_every, folder)
+    istr = get_iter_str(0, t.ncells + 1)
+    base_dir = base_path(folder)
+    sims_dir = sims_path(folder)
+    base_str = string(sims_dir, "output_L=", Int64(grid.domain[2]), "_nc=",
+                      grid.ncells, "_v=", v, "_")
+    dataset = string(base_str, istr, ".h5")
+
+    if !isdir(base_dir)
+        mkdir(base_dir)
+    end
+    if !isdir(sims_dir)
+        mkdir(sims_dir)
+    end
+
+    if isfile(dataset)
+        rm(dataset)
+    end
+
+    println("Output data at directory: ", sims_dir)
+
+    # write initial data
+    xcoord = params.grid_coords
+    pvd = InputOutput.writevtk_initialdata(statevector, sims_dir, xcoord)
+
+
+    println("=========Starting time integration=========")
+    print("Allocating registers for time integrator...")
+    reg2 = similar(statevector)
+    reg3 = similar(statevector)
+    reg4 = similar(statevector)
+    println("✅")
+
+    println("saving every=", save_every)
+
+    nt = t.ncells + 1
+    dt = params.dt
+    indices = CartesianIndices(statevector)
+
+    InputOutput.write_metadata(base_dir, params, save_every, sims_dir)
+    for (i, ti) in enumerate(params.ti)
+        i -= 1
+        if i == 0
+            continue
+        end
+        println("Iteration = ", i, "/", t.ncells)
+        @time RK4(rhs, dt, statevector, reg2, reg3, reg4, params, ti, indices)
+        if i % save_every == 0
+            InputOutput.writevtk(statevector, sims_dir, xcoord, ti, i, pvd)
+        end
+    end
+
+    InputOutput.save_pvd(pvd)
+    return nothing
+end
+
+
+end # end of module