diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 19024d0c..82423681 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,5 +1,5 @@ style = "yas" -margin = 140 +margin = 165 align_assignment = true align_conditional = true whitespace_ops_in_indices = false diff --git a/.buildkite/run_tests.yml b/.buildkite/run_tests.yml index 5339ffb4..19e06b3c 100644 --- a/.buildkite/run_tests.yml +++ b/.buildkite/run_tests.yml @@ -1,3 +1,4 @@ +# DEBUG: until registered, Chmy.jl needs to be explicitly added in "command:" steps steps: - label: "CUDA Julia {{matrix.version}}" matrix: @@ -13,6 +14,7 @@ steps: command: | julia -e 'println("--- :julia: Instantiating project") using Pkg + Pkg.add(url="https://github.com/PTsolvers/Chmy.jl") Pkg.develop(; path=pwd())' || exit 3 julia -e 'println("+++ :julia: Running tests") @@ -38,6 +40,7 @@ steps: command: | julia -e 'println("--- :julia: Instantiating project") using Pkg + Pkg.add(url="https://github.com/PTsolvers/Chmy.jl") Pkg.develop(; path=pwd())' || exit 3 julia -e 'println("+++ :julia: Running tests") @@ -46,7 +49,7 @@ steps: agents: queue: "juliagpu" rocm: "*" - rocmgpu: "gfx1101" # select Ludovic's Navi 3 card (ROCm 6.0) + rocmgpu: "*" #"gfx1101" # select Ludovic's Navi 3 card (ROCm 6.0) timeout_in_minutes: 120 soft_fail: - exit_status: 3 diff --git a/.github/workflows/DocPreviewCleanup.yml b/.github/workflows/DocPreviewCleanup.yml new file mode 100644 index 00000000..2d1d6d18 --- /dev/null +++ b/.github/workflows/DocPreviewCleanup.yml @@ -0,0 +1,34 @@ +name: Doc Preview Cleanup + +on: + pull_request: + types: [closed] + +jobs: + doc-preview-cleanup: + # Do not run on forks to avoid authorization errors + # Source: https://github.community/t/have-github-action-only-run-on-master-repo-and-not-on-forks/140840/18 + # Note: This does not always work as intended - but you can just ignore + # the failed CI runs after merging a PR + if: github.repository_owner == 'PTsolvers' + runs-on: ubuntu-latest + steps: + - name: Checkout gh-pages branch + uses: actions/checkout@v4 + with: + ref: gh-pages + + - name: Delete preview and history + shell: bash + run: | + git config user.name "Documenter.jl" + git config user.email "documenter@juliadocs.github.io" + git rm -rf --ignore-unmatch "previews/PR$PRNUM" + git commit -m "delete preview" --allow-empty + git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree}) + env: + PRNUM: ${{ github.event.number }} + + - name: Push changes + run: | + git push --force origin gh-pages-new:gh-pages diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 31fff6f0..ab8308a6 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -1,3 +1,4 @@ +# DEBUG: until registered, Chmy.jl needs to be explicitly added in "Install dependencies" step name: Documentation on: @@ -14,11 +15,15 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: - version: '1.10' + version: '1' - name: Install dependencies - run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + run: | + julia --project=docs/ -e 'using Pkg + Pkg.add(url="https://github.com/PTsolvers/Chmy.jl") + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate()' - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token diff --git a/.github/workflows/UnitTests.yml b/.github/workflows/UnitTests.yml index 0cb26030..c3c208fb 100644 --- a/.github/workflows/UnitTests.yml +++ b/.github/workflows/UnitTests.yml @@ -1,3 +1,4 @@ +# DEBUG: until registered, Chmy.jl needs to be explicitly added in "Run tests" (instead of using "legacy testing") name: Unit Tests on: push: @@ -45,7 +46,7 @@ jobs: # allow_failure: true steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v1 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} @@ -59,8 +60,21 @@ jobs: ${{ runner.os }}-test-${{ env.cache-name }}- ${{ runner.os }}-test- ${{ runner.os }}- - - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-runtest@latest + # --- legacy testing + # - uses: julia-actions/julia-buildpkg@latest + # - uses: julia-actions/julia-runtest@latest + # --- testing with un-registered Chmy.jl + - name: Run tests + run: | + julia -e 'println("--- :julia: Instantiating project") + using Pkg + Pkg.add(url="https://github.com/PTsolvers/Chmy.jl") + Pkg.develop(; path=pwd())' + + julia -e 'println("+++ :julia: Running tests") + using Pkg + Pkg.test("FastIce"; coverage=true)' + # codecov - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v4 with: diff --git a/Project.toml b/Project.toml index 6d5a66d8..e426f7db 100644 --- a/Project.toml +++ b/Project.toml @@ -1,41 +1,26 @@ name = "FastIce" uuid = "e0de9f13-a007-490e-b696-b07d031015ca" authors = ["Ludovic Raess , Ivan Utkin , Samuel Omlin and contributors"] -version = "0.2.0" +version = "0.3.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -ImplicitGlobalGrid = "4d7a3746-15be-11ea-1130-334b0c4f5fa0" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" -OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -Preferences = "21216c6a-2e73-6563-6e65-726566657250" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -[weakdeps] -AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" - -[extensions] -FastIceAMDGPUExt = "AMDGPU" -FastIceCUDAExt = "CUDA" - [compat] -AMDGPU = "0.8" Adapt = "4" -CUDA = "5" ElasticArrays = "1" GeometryBasics = "0.4" -HDF5 = "0.17" -ImplicitGlobalGrid = "0.15" +HDF5 = "0.16, 0.17" KernelAbstractions = "0.9" LightXML = "0.9" MPI = "0.20" -OffsetArrays = "1" -Preferences = "1" diff --git a/docs/src/lib/modules.md b/docs/src/lib/modules.md index 71d52907..5ef87031 100644 --- a/docs/src/lib/modules.md +++ b/docs/src/lib/modules.md @@ -1,26 +1,5 @@ # Modules -## Grids - -```@autodocs -Modules = [FastIce.Grids] -Order = [:type, :function] -``` - -## Distributed - -```@autodocs -Modules = [FastIce.Distributed] -Order = [:type, :function] -``` - -## Kernel launch - -```@autodocs -Modules = [FastIce.KernelLaunch] -Order = [:type, :function] -``` - ## Writers ```@autodocs @@ -33,4 +12,4 @@ Order = [:type, :function] ```@autodocs Modules = [FastIce.LevelSets] Order = [:type, :function] -``` \ No newline at end of file +``` diff --git a/ext/FastIceAMDGPUExt/FastIceAMDGPUExt.jl b/ext/FastIceAMDGPUExt/FastIceAMDGPUExt.jl deleted file mode 100644 index 30e5a159..00000000 --- a/ext/FastIceAMDGPUExt/FastIceAMDGPUExt.jl +++ /dev/null @@ -1,14 +0,0 @@ -module FastIceAMDGPUExt - -using FastIce, AMDGPU, AMDGPU.ROCKernels -import FastIce.Architectures: heuristic_groupsize, set_device!, get_device - -set_device!(dev::HIPDevice) = AMDGPU.device!(dev) - -get_device(::ROCBackend, id::Integer) = HIPDevice(id) - -heuristic_groupsize(::HIPDevice, ::Val{1}) = (256, ) -heuristic_groupsize(::HIPDevice, ::Val{2}) = (128, 2, ) -heuristic_groupsize(::HIPDevice, ::Val{3}) = (128, 2, 1, ) - -end diff --git a/ext/FastIceCUDAExt/FastIceCUDAExt.jl b/ext/FastIceCUDAExt/FastIceCUDAExt.jl deleted file mode 100644 index 9918672e..00000000 --- a/ext/FastIceCUDAExt/FastIceCUDAExt.jl +++ /dev/null @@ -1,15 +0,0 @@ -module FastIceCUDAExt - -using FastIce, CUDA, CUDA.CUDAKernels - -import FastIce.Architectures: heuristic_groupsize, set_device!, get_device - -set_device!(dev::CuDevice) = CUDA.device!(dev) - -get_device(::CUDABackend, id::Integer) = CuDevice(id - 1) - -heuristic_groupsize(::CuDevice, ::Val{1}) = (256,) -heuristic_groupsize(::CuDevice, ::Val{2}) = (32, 8) -heuristic_groupsize(::CuDevice, ::Val{3}) = (32, 8, 1) - -end diff --git a/scripts_future_API/Project.toml b/scripts_future_API/Project.toml index 9b70abea..94563595 100644 --- a/scripts_future_API/Project.toml +++ b/scripts_future_API/Project.toml @@ -1,9 +1,6 @@ [deps] -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" FastIce = "e0de9f13-a007-490e-b696-b07d031015ca" -FastIceTools = "14532042-4700-46e0-8dec-55d517bc1ae0" -FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" diff --git a/scripts_future_API/bench3d.jl b/scripts_future_API/bench3d.jl deleted file mode 100644 index 3162ca63..00000000 --- a/scripts_future_API/bench3d.jl +++ /dev/null @@ -1,139 +0,0 @@ -using KernelAbstractions -using MPI -using Printf - -using AMDGPU - -# using CUDA -# using NVTX - -include("mpi_utils.jl") -include("mpi_utils2.jl") - -macro inn(A) esc(:( $A[ix+1, iy+1, iz+1] )) end -macro d2_xi(A) esc(:( $A[ix+2, iy+1, iz+1] - $A[ix+1, iy+1, iz+1] - $A[ix+1, iy+1, iz+1] - $A[ix, iy+1, iz+1] )) end -macro d2_yi(A) esc(:( $A[ix+1, iy+2, iz+1] - $A[ix+1, iy+1, iz+1] - $A[ix+1, iy+1, iz+1] - $A[ix+1, iy, iz+1] )) end -macro d2_zi(A) esc(:( $A[ix+1, iy+1, iz+2] - $A[ix+1, iy+1, iz+1] - $A[ix+1, iy+1, iz+1] - $A[ix+1, iy+1, iz] )) end - -@kernel function diffusion_kernel!(A_new, A, h, _dx, _dy, _dz, offset) - ix, iy, iz = @index(Global, NTuple) - ix += offset[1] - 1 - iy += offset[2] - 1 - iz += offset[3] - 1 - if ix ∈ axes(A_new, 1)[2:end-1] && iy ∈ axes(A_new, 2)[2:end-1] && iz ∈ axes(A_new, 3)[2:end-1] - @inbounds A_new[ix, iy, iz] = A[ix, iy, iz] + h * ((A[ix-1, iy , iz ] + A[ix+1, iy , iz ] - 2.0 * A[ix, iy, iz]) * _dx * _dx + - (A[ix , iy-1, iz ] + A[ix , iy+1, iz ] - 2.0 * A[ix, iy, iz]) * _dy * _dy + - (A[ix , iy , iz-1] + A[ix , iy , iz+1] - 2.0 * A[ix, iy, iz]) * _dz * _dz ) - end - # if (ix < size(A, 1) - 2 && iy < size(A, 2) - 2 && iz < size(A, 3) - 2) - # # @inbounds @inn(A_new) = @inn(A) + h - # @inbounds @inn(A_new) = @inn(A) + h * (_dx * _dx * @d2_xi(A) + _dy * _dy * @d2_yi(A) + _dz * _dz * @d2_zi(A)) - # end -end - -function compute_ka(hide_comm, comm, backend, neighbors, ranges, A_new, A, h, _dx, _dy, _dz, iters, me) - (me==0) && print("Starting the time loop 🚀...") - MPI.Barrier(comm) - tic = time_ns() - for _ = 1:iters - # copyto!(A, A_new) - # KernelAbstractions.synchronize(backend) - hide_comm(diffusion_kernel!(backend, 256), neighbors, ranges, A_new, A, h, _dx, _dy, _dz) - A, A_new = A_new, A - - # diffusion_kernel!(backend, 256)(A_new, A, h, _dx, _dy, _dz, (1, 1, 1); ndrange=size(A)) - # KernelAbstractions.synchronize(backend) - # A, A_new = A_new, A - end - wtime = (time_ns() - tic) * 1e-9 - (me==0) && println("done") - return wtime -end - -function main(backend=CPU(), T::DataType=Float64, dims=(0, 0, 0)) - # physics - l = 10.0 - # numerics - iters, warmup = 35, 5 - nx, ny, nz = 1024, 1024, 1024 - b_width = (128, 8, 4) - dims, comm, me, neighbors, coords, device = init_distributed(dims; init_MPI=true) - dx, dy, dz = l ./ (nx, ny, nz) - _dx, _dy, _dz = 1.0 ./ (dx, dy, dz) - h = min(dx, dy ,dz)^2 / 6.1 - # init arrays - x_g = (ix, dx) -> (coords[1] * (nx - 2) + (ix-1)) * dx + dx/2 - y_g = (iy, dy) -> (coords[2] * (ny - 2) + (iy-1)) * dy + dy/2 - z_g = (iz, dz) -> (coords[3] * (nz - 2) + (iz-1)) * dz + dz/2 - # Gaussian - # A_ini = zeros(T, nx, ny, nz) - # A_ini .= [exp(-(x_g(ix, dx) - l/2)^2 - (y_g(iy, dy) - l/2)^2 - (z_g(iz, dz) - l/2)^2) for ix=1:size(A_ini, 1), iy=1:size(A_ini, 2), iz=1:size(A_ini, 3)] - A_ini = rand(T, nx, ny, nz) - - A = KernelAbstractions.allocate(backend, T, nx, ny, nz) - A_new = KernelAbstractions.allocate(backend, T, nx, ny, nz) - KernelAbstractions.copyto!(backend, A, A_ini) - KernelAbstractions.synchronize(backend) - A_new = copy(A) - - ### to be hidden later - ranges = split_ndrange(A, b_width) - - exchangers = ntuple(Val(length(neighbors))) do _ - ntuple(_ -> Exchanger(backend, device), Val(2)) - end - - function hide_comm(f, neighbors, ranges, args...) - f(args..., first(ranges[end]); ndrange=size(ranges[end])) - for dim in reverse(eachindex(neighbors)) - ntuple(Val(2)) do side - rank = neighbors[dim][side] - halo = get_recv_view(Val(side), Val(dim), A_new) - border = get_send_view(Val(side), Val(dim), A_new) - range = ranges[2*(dim-1) + side] - offset, ndrange = first(range), size(range) - start_exchange(exchangers[dim][side], comm, rank, halo, border) do compute_bc - f(args..., offset; ndrange) - if compute_bc - # apply_bcs!(Val(dim), fields, bcs.velocity) - end - KernelAbstractions.synchronize(backend) - end - end - wait.(exchangers[dim]) - end - KernelAbstractions.synchronize(backend) - end - ### to be hidden later - - # GC.gc() - # GC.enable(false) - - compute_ka(hide_comm, comm, backend, neighbors, ranges, A_new, A, h, _dx, _dy, _dz, warmup, me) - - for _ in 1:10 - wtime = compute_ka(hide_comm, comm, backend, neighbors, ranges, A_new, A, h, _dx, _dy, _dz, (iters - warmup), me) - - # GC.enable(true) - # GC.gc() - - MPI.Barrier(comm) - wtime_min = MPI.Allreduce(wtime, MPI.MIN, comm) - wtime_max = MPI.Allreduce(wtime, MPI.MAX, comm) - # perf - A_eff = 2 / 2^30 * (nx-2) * (ny-2) * (nz-2) * sizeof(Float64) - wtime_it = (wtime_min, wtime_max) ./ (iters - warmup) - T_eff = A_eff ./ wtime_it - (me==0) && @printf("Executed %d steps in = %1.3e sec @ T_eff = %1.2f GB/s (max %1.2f) \n", (iters - warmup), wtime_max, round(T_eff[2], sigdigits=3), round(T_eff[1], sigdigits=3)) - # @printf("Executed %d steps in = %1.3e sec (@ T_eff = %1.2f GB/s - device %s) \n", (iters - warmup), wtime, round(T_eff, sigdigits=3), AMDGPU.device_id(AMDGPU.device())) - end - finalize_distributed(; finalize_MPI=true) - return -end - -backend = ROCBackend() -T::DataType = Float64 -dims = (0, 0, 0) - -main(backend, T, dims) -# main() \ No newline at end of file diff --git a/scripts_future_API/generate_synthetic.jl b/scripts_future_API/generate_synthetic.jl new file mode 100644 index 00000000..c18b7427 --- /dev/null +++ b/scripts_future_API/generate_synthetic.jl @@ -0,0 +1,29 @@ +using FastIce.Geometries + +nx, ny = 128, 128 +lx, ly = 5.0, 5.0 +zmin, zmax = 0.0, 1.0 +synth_dem = make_synthetic(nx, ny, lx, ly, zmin, zmax, :turtle) + +@views function visme(synth_dem) + x = LinRange(synth_dem.domain.xmin, synth_dem.domain.xmax, nx + 1) + y = LinRange(synth_dem.domain.ymin, synth_dem.domain.ymax, ny + 1) + bed = synth_dem.z_bed + surf = synth_dem.z_surf + + surf2 = copy(surf) + surf2[surf. Array), + p2_ = plot!(axs.ax2, x_g, dem_ice[:, sly] |> Array), + p3 = heatmap!(axs.ax3, wt_na_c[:, sly, :]; colormap=:turbo), + p4 = heatmap!(axs.ax4, wt_ns_c[:, sly, :]; colormap=:turbo), + p5 = heatmap!(axs.ax5, wt_na_c[slx, :, :]; colormap=:turbo), + p6 = heatmap!(axs.ax6, wt_ns_c[slx, :, :]; colormap=:turbo)) + Colorbar(fig[1, 1][1, 2], plt.p1) + Colorbar(fig[2, 2][1, 2], plt.p4) + # display(fig) + save("levset_$(dims_g[1])_v.png", fig) + end + + if do_h5_save + h5names = String[] + fields = Dict("wt_na" => wt.na.ccc, "wt_ns" => wt.ns.ccc) + outdir = "out_visu_mpi_v" + (me == 0) && mkpath(outdir) + + out_h5 = "results.h5" + (me == 0) && @info "saving HDF5 file" + write_h5(arch, grid, joinpath(outdir, out_h5), fields) + push!(h5names, out_h5) + + (me == 0) && @info "saving XDMF file" + (me == 0) && write_xdmf(arch, grid, joinpath(outdir, "results.xdmf3"), fields, h5names) + end + + return +end + +# main_synthetic(backend) +main_vavilov(backend) + +# MPI.Finalize() diff --git a/scripts_future_API/test_vavilov.jl b/scripts_future_API/test_vavilov.jl new file mode 100644 index 00000000..d4ffa482 --- /dev/null +++ b/scripts_future_API/test_vavilov.jl @@ -0,0 +1,255 @@ +using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.GridOperators, Chmy.Fields, Chmy.KernelLaunch, Chmy.BoundaryConditions +using KernelAbstractions + +using FastIce +using FastIce.Geometries +using FastIce.LevelSets +using FastIce.Physics +using FastIce.Writers + +using FastIce.Models.ImmersedBoundaryFullStokes.Isothermal + +using CUDA +using CairoMakie +using JLD2 + +using Printf +using LinearAlgebra + +using Chmy.Distributed +using MPI +MPI.Init() + +backend = CUDABackend() +do_visu = true +do_h5_save = true + +conv(nx, tx) = tx * ((nx + tx ÷ 2 -1 ) ÷ tx) + +function make_synthetic(arch::Architecture, nx, ny, lx, ly, lz, amp, ω, tanβ, el, gl) + arch_2d = SingleDeviceArchitecture(arch) + grid_2D = UniformGrid(arch_2d; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=(nx, ny)) + + # type = :turtle + generate_ice(x, y, gl, lx, ly) = gl * (1.0 - ((1.7 * x / lx + 0.22)^2 + (0.5 * y / ly)^2)) + generate_bed(x, y, amp, ω, tanβ, el, lx, ly, lz) = lz * (amp * sin(ω * x / lx) * sin(ω * y /ly) + tanβ * x / lx + el + (1.5 * y / ly)^2) + + ice = FunctionField(generate_ice, grid_2D, Vertex(); parameters=(; gl, lx, ly)) + bed = FunctionField(generate_bed, grid_2D, Vertex(); parameters=(; amp, ω, tanβ, el, lx, ly, lz)) + + return (; arch, bed, ice, grid_2D, lx, ly, lz) +end + +function main_synthetic(backend=CPU()) + # set-up distributed + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + + # synthetic topo + lx, ly, lz = 5.0, 5.0, 1.0 + amp = 0.1 + ω = 10π + tanβ = tan(-π/6) + el = 0.35 + gl = 0.9 + + nx, ny = 126, 126 + nz = max(conv(ceil(Int, lz / lx * nx), 30), 30) + resol = (nx, ny, nz) + + synthetic_elevation = make_synthetic(arch, nx, ny, lx, ly, lz, amp, ω, tanβ, el, gl) + + run_simulation(synthetic_elevation..., resol...) + return +end + +function extract_dem(arch::Architecture, data_path::String) + data = load(data_path) + dtype = first(keys(data)) + + z_ice = data[dtype].z_surf + z_bed = data[dtype].z_bed + dm = data[dtype].domain + dm = dilate(dm, 0.05) + lx, ly, lz = extents(dm) + nx, ny = size(z_ice) .- 1 + + z_ice .-= dm.zmin # put Z-origin at 0 + z_bed .-= dm.zmin # put Z-origin at 0 + + # rescale tmp + z_ice ./= lz + z_bed ./= lz + lx /= lx / 10 + ly /= ly / 10 + lz /= lz + + arch_2d = SingleDeviceArchitecture(arch) + grid_2D = UniformGrid(arch_2d; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=(nx, ny)) + + ice = Field(arch_2d, grid_2D, Vertex()) + bed = Field(arch_2d, grid_2D, Vertex()) + set!(ice, z_ice) + set!(bed, z_bed) + + return (; arch, bed, ice, grid_2D, lx, ly, lz) +end + +function main_vavilov(backend=CPU()) + # load dem data + data_path = "./vavilov_dem2.jld2" + + # set-up distributed + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + + data_elevation = extract_dem(arch, data_path) + + # nx, ny = 126, 126 + # nz = max(conv(ceil(Int, data_elevation.lz / data_elevation.lx * (nx + 2)), 32), 62) - 2 + # resol = (nx, ny, nz) + resol = (510, 510, 126) + + run_simulation(data_elevation..., resol...) + return +end + +# function main(backend; res) +function run_simulation(arch::Architecture, bed, ice, grid_2D, lx, ly, lz, nx, ny, nz) + topo = topology(arch) + me = global_rank(topo) + # geometry + (me == 0) && @info "size = ($nx, $ny, $nz), extent = ($lx, $ly, $lz)" + dims_l = (nx, ny, nz) + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) + launch = Launcher(arch, grid; outer_width=(16, 8, 4)) + # compute level sets + ψ = (ns=Field(arch, grid, Vertex()), + na=Field(arch, grid, Vertex())) + + compute_levelset_from_dem!(arch, launch, ψ.na, ice, grid_2D, grid) + compute_levelset_from_dem!(arch, launch, ψ.ns, bed, grid_2D, grid) + + invert_levelset!(arch, launch, ψ.ns, grid) + + free_slip = BoundaryCondition{Traction}(0.0, 0.0, 0.0) + + boundary_conditions = (x=(free_slip, free_slip), + y=(free_slip, free_slip), + z=(free_slip, free_slip)) + + gravity = (x=ValueField(0.0), + y=ValueField(0.0), + z=ValueField(1.0)) + + # numerics + niter = 50maximum(size(grid, Center())) + ncheck = 5maximum(size(grid, Center())) + + r = 0.9 + re_mech = 5π + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.1) + θ_dτ = lτ_re_m * (r + 4 / 3) / vdτ + dτ_r = 1.0 / (θ_dτ + 1.0) + nudτ = vdτ * lτ_re_m + + solver_params = (Δτ=(Pr=r / θ_dτ, τ=(xx=dτ_r, yy=dτ_r, zz=dτ_r, xy=dτ_r, xz=dτ_r, yz=dτ_r), V=(x=nudτ, y=nudτ, z=nudτ)),) + + rheology = LinearViscousRheology(OneField{Float64}()) + + model = IsothermalImmersedBoundaryFullStokesModel(; + arch, + grid, + boundary_conditions, + gravity, + rheology, + solver_params, + level_sets=ψ) + + # init + ω_from_ψ!(arch, launch, model.field_masks.ns, ψ.ns, grid) + ω_from_ψ!(arch, launch, model.field_masks.na, ψ.na, grid) + + set!(model.viscosity.η, rheology.η) + set!(model.viscosity.η_next, rheology.η) + + bc!(arch, grid, model.viscosity.η => Neumann()) + bc!(arch, grid, model.viscosity.η_next => Neumann()) + + nx, ny, nz = size(grid, Center()) + + iy_sl = ny ÷ 2 + + Vm = Field(arch, grid, Center()) + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + + fig = Figure(; size=(800, 400)) + axs = (Pr = Axis(fig[1, 1][1, 1]; aspect=DataAspect(), ylabel="z", title="Pr"), + Vm = Axis(fig[1, 2][1, 1]; aspect=DataAspect(), title="|V|"), + ωna = Axis(fig[2, 1][1, 1]; aspect=DataAspect(), xlabel="x"), + ωns = Axis(fig[2, 2][1, 1]; aspect=DataAspect(), xlabel="x")) + plt = (Pr = heatmap!(axs.Pr, xcenters(grid), zcenters(grid), Array(interior(model.stress.P)[:, iy_sl, :]); colormap=:turbo), + Vm = heatmap!(axs.Vm, xcenters(grid), zcenters(grid), Array(interior(Vm)[:, iy_sl, :]); colormap=:turbo), + ωna = heatmap!(axs.ωna, xcenters(grid), zcenters(grid), Array(interior(model.field_masks.na.ccc)[:, iy_sl, :]); colormap=:turbo), + ωns = heatmap!(axs.ωns, xcenters(grid), zcenters(grid), Array(interior(model.field_masks.ns.ccc)[:, iy_sl, :]); colormap=:turbo)) + Colorbar(fig[1, 1][1, 2], plt.Pr) + Colorbar(fig[1, 2][1, 2], plt.Vm) + Colorbar(fig[2, 1][1, 2], plt.ωna) + Colorbar(fig[2, 2][1, 2], plt.ωns) + + display(fig) + + if do_h5_save + h5names = String[] + fields = Dict("Pr" => model.stress.P, "Vm" => Vm, "wt_ns_c" => model.field_masks.ns.ccc, "wt_na_c" => model.field_masks.na.ccc) + outdir = "out_visu" + mkpath(outdir) + end + + iter = 1 + for iter in 1:niter + advance_iteration!(model, 0.0, 1.0) + if (iter % ncheck == 0) + compute_residuals!(model) + err = (Pr = norm(model.residual.r_P, Inf), + Vx = norm(model.residual.r_V.x, Inf), + Vy = norm(model.residual.r_V.y, Inf), + Vz = norm(model.residual.r_V.z, Inf)) + if any(.!isfinite.(values(err))) + error("simulation failed, err = $err") + end + iter_nx = iter / max(nx, ny, nz) + @printf(" iter/nx = %.1f, err = [Pr = %1.3e, Vx = %1.3e, Vy = %1.3e, Vz = %1.3e]\n", iter_nx, err...) + end + end + if do_visu + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + + plt.Pr[3][] = interior(model.stress.P)[:, iy_sl, :] + plt.Vm[3][] = interior(Vm)[:, iy_sl, :] + yield() + display(fig) + end + if do_h5_save + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + out_h5 = "results.h5" + @info "saving HDF5 file" + write_h5(arch, grid, joinpath(outdir, out_h5), fields) + push!(h5names, out_h5) + + @info "saving XDMF file" + write_xdmf(arch, grid, joinpath(outdir, "results.xdmf3"), fields, h5names) + end + return +end + +# main_synthetic(backend) +main_vavilov(backend) + +# MPI.Finalize() diff --git a/scripts_future_API/test_volume_fractions.jl b/scripts_future_API/test_volume_fractions.jl new file mode 100644 index 00000000..ca0f9e9d --- /dev/null +++ b/scripts_future_API/test_volume_fractions.jl @@ -0,0 +1,175 @@ +using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.GridOperators, Chmy.Fields, Chmy.KernelLaunch, Chmy.BoundaryConditions +using KernelAbstractions + +using FastIce +using FastIce.LevelSets +using FastIce.Physics +using FastIce.Models.ImmersedBoundaryFullStokes.Isothermal + +using FastIce.Writers + +# using CUDA +# using GLMakie +using CairoMakie + +using Printf +using LinearAlgebra + +function main(backend; res) + arch = Arch(backend) + grid = UniformGrid(arch; origin=(-1, -1, 0), extent=(2, 2, 1), dims=res) + + grid_2D = UniformGrid(arch; origin=(-1, -1), extent=(2, 2), dims=(100, 100)) + + # bed parameters + amp = 0.05 + ωb = 10π / 2.0 + αx = -tan(π / 20) + αy = -tan(π / 25) + z0_bed = 0.3 + + # ice parameters + x0 = -0.1 + y0 = -0.1 + z0 = -0.1 + rad = 0.9 + + bed = FunctionField(grid_2D, Vertex(); parameters=(amp, ωb, αx, αy, z0_bed)) do x, y, amp, ωb, αx, αy, z0_bed + return z0_bed + x * αx + y * αy + amp * sin(ωb * x) * cos(ωb * y) + end + + ice = FunctionField(grid_2D, Vertex(); parameters=(x0, y0, z0, rad)) do x, y, x0, y0, z0, rad + return z0 + sqrt(max(rad^2 - (x - x0)^2 - (y - y0)^2, 0.0)) + end + + ψ = (ns=Field(arch, grid, Vertex()), + na=Field(arch, grid, Vertex())) + + launch = Launcher(arch, grid) + + compute_levelset_from_dem!(arch, launch, ψ.na, ice, grid_2D, grid) + compute_levelset_from_dem!(arch, launch, ψ.ns, bed, grid_2D, grid) + + invert_levelset!(arch, launch, ψ.ns, grid) + + free_slip = BoundaryCondition{Traction}(0.0, 0.0, 0.0) + + boundary_conditions = (x=(free_slip, free_slip), + y=(free_slip, free_slip), + z=(free_slip, free_slip)) + + gravity = (x=ValueField(0.0), + y=ValueField(0.0), + z=ValueField(1.0)) + + # numerics + niter = 25maximum(size(grid, Center())) + ncheck = 5maximum(size(grid, Center())) + do_visu = true + do_h5_save = true + + r = 0.9 + re_mech = 5π + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.1) + θ_dτ = lτ_re_m * (r + 4 / 3) / vdτ + dτ_r = 1.0 / (θ_dτ + 1.0) + nudτ = vdτ * lτ_re_m + + solver_params = (Δτ=(Pr=r / θ_dτ, τ=(xx=dτ_r, yy=dτ_r, zz=dτ_r, xy=dτ_r, xz=dτ_r, yz=dτ_r), V=(x=nudτ, y=nudτ, z=nudτ)),) + + rheology = LinearViscousRheology(OneField{Float64}()) + + model = IsothermalImmersedBoundaryFullStokesModel(; + arch, + grid, + boundary_conditions, + gravity, + rheology, + solver_params, + level_sets=ψ) + + # init + ω_from_ψ!(arch, launch, model.field_masks.ns, ψ.ns, grid) + ω_from_ψ!(arch, launch, model.field_masks.na, ψ.na, grid) + + set!(model.viscosity.η, rheology.η) + set!(model.viscosity.η_next, rheology.η) + + bc!(arch, grid, model.viscosity.η => Neumann()) + bc!(arch, grid, model.viscosity.η_next => Neumann()) + + nx, ny, nz = size(grid, Center()) + + iy_sl = ny ÷ 2 + + Vm = Field(arch, grid, Center()) + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + + fig = Figure(; size=(800, 400)) + axs = (Pr = Axis(fig[1, 1][1, 1]; aspect=DataAspect(), ylabel="z", title="Pr"), + Vm = Axis(fig[1, 2][1, 1]; aspect=DataAspect(), title="|V|"), + ωna = Axis(fig[2, 1][1, 1]; aspect=DataAspect(), xlabel="x"), + ωns = Axis(fig[2, 2][1, 1]; aspect=DataAspect(), xlabel="x")) + plt = (Pr = heatmap!(axs.Pr, xcenters(grid), zcenters(grid), Array(interior(model.stress.P)[:, iy_sl, :]); colormap=:turbo), + Vm = heatmap!(axs.Vm, xcenters(grid), zcenters(grid), Array(interior(Vm)[:, iy_sl, :]); colormap=:turbo), + ωna = heatmap!(axs.ωna, xcenters(grid), zcenters(grid), Array(interior(model.field_masks.na.ccc)[:, iy_sl, :]); colormap=:turbo), + ωns = heatmap!(axs.ωns, xcenters(grid), zcenters(grid), Array(interior(model.field_masks.ns.ccc)[:, iy_sl, :]); colormap=:turbo)) + Colorbar(fig[1, 1][1, 2], plt.Pr) + Colorbar(fig[1, 2][1, 2], plt.Vm) + Colorbar(fig[2, 1][1, 2], plt.ωna) + Colorbar(fig[2, 2][1, 2], plt.ωns) + + display(fig) + + if do_h5_save + h5names = String[] + fields = Dict("Pr" => model.stress.P, "Vm" => Vm, "wt_ns_c" => model.field_masks.ns.ccc, "wt_na_c" => model.field_masks.na.ccc) + outdir = "out_visu" + mkpath(outdir) + end + + iter = 1 + for iter in 1:niter + advance_iteration!(model, 0.0, 1.0) + if (iter % ncheck == 0) + compute_residuals!(model) + err = (Pr = norm(model.residual.r_P, Inf), + Vx = norm(model.residual.r_V.x, Inf), + Vy = norm(model.residual.r_V.y, Inf), + Vz = norm(model.residual.r_V.z, Inf)) + if any(.!isfinite.(values(err))) + error("simulation failed, err = $err") + end + iter_nx = iter / max(nx, ny, nz) + @printf(" iter/nx = %.1f, err = [Pr = %1.3e, Vx = %1.3e, Vy = %1.3e, Vz = %1.3e]\n", iter_nx, err...) + end + end + if do_visu + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + + plt.Pr[3][] = interior(model.stress.P)[:, iy_sl, :] + plt.Vm[3][] = interior(Vm)[:, iy_sl, :] + yield() + display(fig) + end + if do_h5_save + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + out_h5 = "results.h5" + @info "saving HDF5 file" + write_h5(arch, grid, joinpath(outdir, out_h5), fields) + push!(h5names, out_h5) + + @info "saving XDMF file" + write_xdmf(arch, grid, joinpath(outdir, "results.xdmf3"), fields, h5names) + end + return +end + +main(CPU(); res=(128, 128, 64) .- 2) diff --git a/scripts_future_API/test_volume_fractions_mpi.jl b/scripts_future_API/test_volume_fractions_mpi.jl new file mode 100644 index 00000000..4323cb8e --- /dev/null +++ b/scripts_future_API/test_volume_fractions_mpi.jl @@ -0,0 +1,199 @@ +using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.GridOperators, Chmy.Fields, Chmy.KernelLaunch, Chmy.BoundaryConditions +using KernelAbstractions + +using FastIce +using FastIce.LevelSets +using FastIce.Physics +using FastIce.Models.ImmersedBoundaryFullStokes.Isothermal + +using FastIce.Writers + +# using GLMakie +using CairoMakie + +backend = CPU() + +# using AMDGPU +# backend = ROCBackend() + +using Printf +using LinearAlgebra + +using Chmy.Distributed +using MPI +MPI.Init() + +function max_mpi(A) + max_l = maximum(A) + return MPI.Allreduce(max_l, MPI.MAX, MPI.COMM_WORLD) +end + +function main(backend=CPU(); res) + + # 3D distributed§ + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + topo = topology(arch) + me = global_rank(topo) + dims_l = res + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-1, -1, 0), extent=(2, 2, 1), dims=dims_g) + + # 2D single device + arch_2D = SingleDeviceArchitecture(arch) + grid_2D = UniformGrid(arch_2D; origin=(-1, -1), extent=(2, 2), dims=(100, 100)) + + # bed parameters + amp = 0.05 + ωb = 10π / 2.0 + αx = -tan(π / 20) + αy = -tan(π / 25) + z0_bed = 0.3 + + # ice parameters + x0 = -0.1 + y0 = -0.1 + z0 = -0.1 + rad = 0.9 + + bed = FunctionField(grid_2D, Vertex(); parameters=(amp, ωb, αx, αy, z0_bed)) do x, y, amp, ωb, αx, αy, z0_bed + return z0_bed + x * αx + y * αy + amp * sin(ωb * x) * cos(ωb * y) + end + + ice = FunctionField(grid_2D, Vertex(); parameters=(x0, y0, z0, rad)) do x, y, x0, y0, z0, rad + return z0 + sqrt(max(rad^2 - (x - x0)^2 - (y - y0)^2, 0.0)) + end + + ψ = (ns=Field(arch, grid, Vertex()), + na=Field(arch, grid, Vertex())) + + free_slip = BoundaryCondition{Slip}(0.0, 0.0, 0.0) + + boundary_conditions = (x=(free_slip, free_slip), + y=(free_slip, free_slip), + z=(free_slip, free_slip)) + + gravity = (x=ValueField(0.0), + y=ValueField(0.0), + z=ValueField(1.0)) + + # numerics + niter = 25maximum(size(grid, Center())) + ncheck = 5maximum(size(grid, Center())) + do_visu = true + do_h5_save = false + + r = 0.9 + re_mech = 5π + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.1) + θ_dτ = lτ_re_m * (r + 4 / 3) / vdτ + dτ_r = 1.0 / (θ_dτ + 1.0) + nudτ = vdτ * lτ_re_m + + solver_params = (Δτ=(Pr=r / θ_dτ, τ=(xx=dτ_r, yy=dτ_r, zz=dτ_r, xy=dτ_r, xz=dτ_r, yz=dτ_r), V=(x=nudτ, y=nudτ, z=nudτ)),) + + rheology = LinearViscousRheology(OneField{Float64}()) + + model = IsothermalImmersedBoundaryFullStokesModel(; + arch, + grid, + boundary_conditions, + gravity, + rheology, + solver_params, + level_sets=ψ) + + # init + compute_levelset_from_dem!(arch, model.launcher, ψ.na, ice, grid_2D, grid) + compute_levelset_from_dem!(arch, model.launcher, ψ.ns, bed, grid_2D, grid) + + invert_levelset!(arch, model.launcher, ψ.ns, grid) + + ω_from_ψ!(arch, model.launcher, model.field_masks.ns, ψ.ns, grid) # exchange_halo in ω_from_ψ! wrapper for DistributedArch + ω_from_ψ!(arch, model.launcher, model.field_masks.na, ψ.na, grid) + + set!(model.viscosity.η, rheology.η) + set!(model.viscosity.η_next, rheology.η) + + bc!(arch, grid, model.viscosity.η => Neumann(); exchange=model.viscosity.η) + bc!(arch, grid, model.viscosity.η_next => Neumann(); exchange=model.viscosity.η_next) + + iter = 1 + for iter in 1:niter + advance_iteration!(model, 0.0, 1.0) + if (iter % ncheck == 0) + compute_residuals!(model) + err = (Pr = max_mpi(abs.(interior(model.residual.r_P ))), + Vx = max_mpi(abs.(interior(model.residual.r_V.x))), + Vy = max_mpi(abs.(interior(model.residual.r_V.y))), + Vz = max_mpi(abs.(interior(model.residual.r_V.z)))) + if any(.!isfinite.(values(err))) + error("simulation failed, err = $err") + end + iter_nx = iter / maximum(size(grid, Center())) + (me == 0) && @printf(" iter/nx = %.1f, err = [Pr = %1.3e, Vx = %1.3e, Vy = %1.3e, Vz = %1.3e]\n", iter_nx, err...) + end + end + if do_visu + Vm = Field(arch, grid, Center()) + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + + d1 = model.field_masks.ns.ccc + d2 = model.field_masks.na.ccc + d3 = Vm + d4 = model.stress.P + tp1 = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(d1)) .* dims(topo)) : nothing + tp2 = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(d2)) .* dims(topo)) : nothing + tp3 = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(d3)) .* dims(topo)) : nothing + tp4 = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(d4)) .* dims(topo)) : nothing + gather!(arch, tp1, d1) + gather!(arch, tp2, d2) + gather!(arch, tp3, d3) + gather!(arch, tp4, d4) + if me == 0 + slx = ceil(Int, size(tp1, 1) / 2) + sly = ceil(Int, size(tp1, 2) / 2) + slz = ceil(Int, size(tp1, 3) / 2) + fig = Figure(; size=(1000, 500), fontsize=22) + axs = (ax1 = Axis(fig[1, 1]; aspect=DataAspect(), title="ns.ccc"), + ax2 = Axis(fig[1, 2]; aspect=DataAspect(), title="na.ccc"), + ax3 = Axis(fig[2, 1]; aspect=DataAspect(), title="Vm"), + ax4 = Axis(fig[2, 2]; aspect=DataAspect(), title="Pr")) + plt = (p1 = heatmap!(axs.ax1, tp1[:, sly, :]; colormap=:turbo), + p2 = heatmap!(axs.ax2, tp2[:, sly, :]; colormap=:turbo), + p3 = heatmap!(axs.ax3, tp3[:, sly, :]; colormap=:turbo), + p4 = heatmap!(axs.ax4, tp4[:, sly, :]; colormap=:turbo)) + Colorbar(fig[1, 1][1, 2], plt.p1) + Colorbar(fig[1, 2][1, 2], plt.p2) + Colorbar(fig[2, 1][1, 2], plt.p3) + Colorbar(fig[2, 2][1, 2], plt.p4) + save("out_vis.png", fig) + end + end + if do_h5_save + h5names = String[] + fields = Dict("Pr" => model.stress.P, "Vm" => Vm, "wt_ns_c" => model.field_masks.ns.ccc, "wt_na_c" => model.field_masks.na.ccc) + outdir = "out_visu" + (me == 0) && mkpath(outdir) + + set!(Vm, grid, (g, loc, ix, iy, iz, V) -> sqrt(lerp(V.x, loc, g, ix, iy, iz)^2 + + lerp(V.y, loc, g, ix, iy, iz)^2 + + lerp(V.z, loc, g, ix, iy, iz)^2); discrete=true, parameters=(model.velocity.V,)) + + out_h5 = "results.h5" + (me == 0) && @info "saving HDF5 file" + write_h5(arch, grid, joinpath(outdir, out_h5), fields) + push!(h5names, out_h5) + + (me == 0) && @info "saving XDMF file" + (me == 0) && write_xdmf(arch, grid, joinpath(outdir, "results.xdmf3"), fields, h5names) + end + return +end + +# main(backend; res=(128, 128, 64) .- 2) +main(backend; res=(32, 32, 32) .- 2) + +MPI.Finalize() diff --git a/scripts_future_API/tm_stokes_chmy.jl b/scripts_future_API/tm_stokes_chmy.jl new file mode 100644 index 00000000..e60c56c2 --- /dev/null +++ b/scripts_future_API/tm_stokes_chmy.jl @@ -0,0 +1,118 @@ +using Chmy.Architectures +using Chmy.Grids +using Chmy.Fields +using Chmy.BoundaryConditions + +using FastIce, FastIce.Models.FullStokes.Isothermal, FastIce.Physics + +using KernelAbstractions + +using LinearAlgebra +using Printf + +using CairoMakie + +backend = CPU() +arch = Arch(backend) + +outer_width = (4, 4, 4) #(128, 32, 4)# + +grid = UniformGrid(arch; + origin=(-0.5, -0.5, 0.0), + extent=(1.0, 1.0, 1.0), + dims=(62, 62, 62)) + +FastIce.greet_fast(; bold=true, color=:blue) + +const VBC = BoundaryCondition{Velocity} +const SBC = BoundaryCondition{Slip} +const TBC = BoundaryCondition{Traction} + +free_slip = SBC(0.0, 0.0, 0.0) +free_surface = TBC(0.0, 0.0, 0.0) +no_slip = VBC(0.0, 0.0, 0.0) + +boundary_conditions = (x=(free_slip, free_slip), + y=(free_slip, free_slip), + z=(no_slip, free_surface)) + +gravity = (x=ZeroField{Float64}(), + y=ZeroField{Float64}(), + z=ValueField(1.0)) + +# numerics +niter = 100maximum(size(grid, Center())) +ncheck = 2maximum(size(grid, Center())) +do_visu = true + +r = 0.7 +re_mech = 2π +lτ_re_m = minimum(extent(grid, Vertex())) / re_mech +vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.5) +θ_dτ = lτ_re_m * (r + 4 / 3) / vdτ +dτ_r = 1.0 / (θ_dτ + 1.0) +nudτ = vdτ * lτ_re_m + +solver_params = (Δτ=(Pr=r / θ_dτ, τ=(xx=dτ_r, yy=dτ_r, zz=dτ_r, xy=dτ_r, xz=dτ_r, yz=dτ_r), V=(x=nudτ, y=nudτ, z=nudτ)),) + +init_incl(x, y, z, x0, y0, z0, r, ηi, ηm) = ifelse((x - x0)^2 + (y - y0)^2 + (z - z0)^2 < r^2, ηi, ηm) + +η = FunctionField(init_incl, grid, Center(); parameters=(x0=0.0, y0=0.0, z0=0.5, r=0.1, ηi=1e-1, ηm=1.0)) +rheology = LinearViscousRheology(η) + +model = IsothermalFullStokesModel(; + arch, + grid, + boundary_conditions, + gravity, + rheology, + solver_params, + outer_width) + +fig = Figure() +axs = (Pr = Axis(fig[1, 1][1, 1]; aspect=DataAspect(), xlabel="x", ylabel="z", title="Pr"), + Vx = Axis(fig[1, 2][1, 1]; aspect=DataAspect(), xlabel="x", ylabel="z", title="Vx"), + Vy = Axis(fig[2, 1][1, 1]; aspect=DataAspect(), xlabel="x", ylabel="z", title="Vy"), + Vz = Axis(fig[2, 2][1, 1]; aspect=DataAspect(), xlabel="x", ylabel="z", title="Vz")) +plt = (Pr = heatmap!(axs.Pr, xcenters(grid), zcenters(grid), Array(interior(model.stress.P)[:, size(grid, Center(), Val(2))÷2+1, :]); colormap=:turbo), + Vx = heatmap!(axs.Vx, xvertices(grid), zcenters(grid), Array(interior(model.velocity.V.x)[:, size(grid, Center(), Val(2))÷2+1, :]); colormap=:turbo), + Vy = heatmap!(axs.Vy, xcenters(grid), zcenters(grid), Array(interior(model.velocity.V.y)[:, size(grid, Vertex(), Val(2))÷2+1, :]); colormap=:turbo), + Vz = heatmap!(axs.Vz, xcenters(grid), zvertices(grid), Array(interior(model.velocity.V.z)[:, size(grid, Center(), Val(2))÷2+1, :]); colormap=:turbo)) +Colorbar(fig[1, 1][1, 2], plt.Pr) +Colorbar(fig[1, 2][1, 2], plt.Vx) +Colorbar(fig[2, 1][1, 2], plt.Vy) +Colorbar(fig[2, 2][1, 2], plt.Vz) + +set!(model.stress.P, 0.0) +set!(model.stress.τ, 0.0) +set!(model.velocity.V, 0.0) + +set!(model.viscosity.η, η) +set!(model.viscosity.η_next, η) + +bc!(arch, grid, model.viscosity.η => Neumann()) +bc!(arch, grid, model.viscosity.η_next => Neumann()) + +for iter in 1:niter + advance_iteration!(model, 0.0, 1.0) + if (iter % ncheck == 0) + compute_residuals!(model) + err = (Pr = norm(model.residual.r_P, Inf), + Vx = norm(model.residual.r_V.x, Inf), + Vy = norm(model.residual.r_V.y, Inf), + Vz = norm(model.residual.r_V.z, Inf)) + if any(.!isfinite.(values(err))) + error("simulation failed, err = $err") + end + iter_nx = iter / maximum(size(grid, Center())) + @printf(" iter/nx = %.1f, err = [Pr = %1.3e, Vx = %1.3e, Vy = %1.3e, Vz = %1.3e]\n", iter_nx, err...) + if do_visu + plt.Pr[3][] = interior(model.stress.P)[:, size(grid, Center(), Val(2))÷2+1, :] + plt.Vx[3][] = interior(model.velocity.V.x)[:, size(grid, Center(), Val(2))÷2+1, :] + plt.Vy[3][] = interior(model.velocity.V.y)[:, size(grid, Vertex(), Val(2))÷2+0, :] + plt.Vz[3][] = interior(model.velocity.V.z)[:, size(grid, Center(), Val(2))÷2+1, :] + # yield() + display(fig) + end + end +end diff --git a/scripts_future_API/tm_stokes_mpi_wip.jl b/scripts_future_API/tm_stokes_mpi_wip.jl index 2ca6d143..f18d8785 100644 --- a/scripts_future_API/tm_stokes_mpi_wip.jl +++ b/scripts_future_API/tm_stokes_mpi_wip.jl @@ -1,12 +1,6 @@ using FastIce -using FastIce.Architectures -using FastIce.Grids -using FastIce.Fields -using FastIce.Utils -using FastIce.BoundaryConditions using FastIce.Models.FullStokes.Isothermal using FastIce.Physics -using FastIce.KernelLaunch using FastIce.Writers using FastIce.Logging @@ -18,78 +12,58 @@ using LinearAlgebra, Printf using KernelAbstractions # using AMDGPU -using FastIce.Distributed -using MPI using Logging +using MPI +MPI.Init() using CairoMakie norm_g(A) = (sum2_l = sum(interior(A) .^ 2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD))) max_abs_g(A) = (max_l = maximum(abs.(interior(A))); MPI.Allreduce(max_l, MPI.MAX, MPI.COMM_WORLD)) -@views avx(A) = 0.5 .* (A[1:end-1, :, :] .+ A[2:end, :, :]) -@views avy(A) = 0.5 .* (A[:, 1:end-1, :] .+ A[:, 2:end, :]) -@views avz(A) = 0.5 .* (A[:, :, 1:end-1] .+ A[:, :, 2:end]) - -@views av_xy(A) = 0.25 .* (A[1:end-1, 1:end-1, :] .+ A[2:end, 1:end-1, :] .+ A[2:end, 2:end, :] .+ A[1:end-1, 2:end, :]) -@views av_xz(A) = 0.25 .* (A[1:end-1, :, 1:end-1] .+ A[2:end, :, 1:end-1, :] .+ A[2:end, :, 2:end, :] .+ A[1:end-1, :, 2:end]) -@views av_yz(A) = 0.25 .* (A[:, 1:end-1, 1:end-1] .+ A[:, 2:end, 1:end-1] .+ A[:, 2:end, 2:end] .+ A[:, 1:end-1, 2:end]) - -function main(; do_visu=false, do_save=false, do_h5_save=false) - MPI.Init() - - backend = CPU() - # dims = (4, 2, 2) - # dims = (4, 2, 2) - dims = (0, 0, 0) - topo = CartesianTopology(dims) - arch = Architecture(backend, topo) - set_device!(arch) - - comm = cartesian_communicator(topo) - me = global_rank(topo) # rank - - size_l = (14, 14, 14) - size_g = global_grid_size(topo, size_l) - - outer_width = (4, 4, 4) #(128, 32, 4)# - - grid_g = CartesianGrid(; origin=(-2.0, -1.0, 0.0), - extent=(4.0, 2.0, 2.0), - size=size_g) - - grid_l = local_grid(grid_g, topo) +function main(backend=CPU(); dims_l=(126, 126, 126), do_visu=false, do_save=false, do_h5_save=false) + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + topo = topology(arch) + me = global_rank(topo) + # geometry + lx, ly, lz = 2.0, 2.0, 2.0 + # mechanics + η = 1.0e1 # solid viscosity + G = 1.0e0 # shear modulus + # gravity + # gravity force density + ρg = (x=0.0, y=0.0, z=1.0) + # scales + psc = ρg.z * lz + τsc = η / psc + # numerics + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-lx / 2, -ly / 2, -lz / 2), extent=(lx, ly, lz), dims=dims_g) + launch = Launcher(arch, grid; outer_width=(4, 4, 4)) - # global_logger(FastIce.Logging.MPILogger(0, comm, global_logger())) + (me == 0) && FastIce.greet_fast(; bold=true, color=:blue) - if me == 0 - FastIce.greet_fast(bold=true, color=:blue) - printstyled(grid_g; bold=true) - end - - no_slip = VBC(0.0, 0.0, 0.0) - free_slip = SBC(0.0, 0.0, 0.0) - free_surface = TBC(0.0, 0.0, 0.0) + free_slip = SBC(0.0, 0.0, 0.0) boundary_conditions = (x=(free_slip, free_slip), y=(free_slip, free_slip), - z=(no_slip, free_surface)) + z=(free_slip, free_slip)) + + init_incl(x, y, z, x0, y0, z0, r, in, out) = ifelse((x - x0)^2 + (y - y0)^2 + (z - z0)^2 < r^2, in, out) + ρgz = FunctionField(init_incl, grid, (Center(), Center(), Vertex()); parameters=(x0=0.0, y0=0.0, z0=0.0, r=0.1lx, in=ρg.z, out=0.0)) - ρgx(x, y, z) = 0.25 - ρgy(x, y, z) = 0.0 - ρgz(x, y, z) = 1.0 - gravity = (x=FunctionField(ρgx, grid_l, (Vertex(), Center(), Center())), - y=FunctionField(ρgy, grid_l, (Center(), Vertex(), Center())), - z=FunctionField(ρgz, grid_l, (Center(), Center(), Vertex()))) + gravity = (x=ZeroField{Float64}(), + y=ZeroField{Float64}(), + z=FunctionField(ρgz, grid, (Center(), Center(), Vertex()))) # numerics - niter = 20maximum(size(grid_g)) - ncheck = 4maximum(size(grid_g)) + niter = 20maximum(size(grid, Center())) + ncheck = 4maximum(size(grid, Center())) r = 0.7 re_mech = 4π - lτ_re_m = minimum(extent(grid_g)) / re_mech - vdτ = minimum(spacing(grid_g)) / sqrt(ndims(grid_g) * 1.5) + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.5) θ_dτ = lτ_re_m * (r + 4 / 3) / vdτ dτ_r = 1.0 / (θ_dτ + 1.0) nudτ = vdτ * lτ_re_m @@ -218,16 +192,6 @@ function main(; do_visu=false, do_save=false, do_h5_save=false) (me == 0) && write_xdmf(arch, grid_g, joinpath(outdir, "results.xdmf3"), fields, h5names) end - if do_h5_save - out_h5 = "results.h5" - (me == 0) && @info "saving HDF5 file" - write_h5(arch, grid_g, joinpath(outdir, out_h5), fields) - push!(h5names, out_h5) - - (me == 0) && @info "saving XDMF file" - (me == 0) && write_xdmf(arch, grid_g, joinpath(outdir, "results.xdmf3"), fields, h5names) - end - if do_save || do_visu copyto!(Pr_v, interior(model.stress.Pr)) copyto!(τxx_v, interior(model.stress.τ.xx)) diff --git a/scripts_future_API/visme_ms_3d.jl b/scripts_future_API/visme_ms_3d.jl new file mode 100644 index 00000000..f54c1cb2 --- /dev/null +++ b/scripts_future_API/visme_ms_3d.jl @@ -0,0 +1,79 @@ +using GLMakie + +const ω = 1π + +# manufactured solution for the confined Stokes flow with free-surface boundaries +# velocity +# manufactured solution for the confined Stokes flow with free-slip boundaries +# helper functions +f(ξ, η) = cos(π * ξ) * (η^2 - 1)^2 - + cos(π * η) * (ξ^2 - 1)^2 +g(ξ, η) = sin(π * η) * (ξ^2 - 1) * ξ - + sin(π * ξ) * (η^2 - 1) * η +p(ξ, η) = cos(π * η) * (1 - 3 * ξ^2) * 2 - + cos(π * ξ) * (1 - 3 * η^2) * 2 +# velocity +vx(x, y, z) = sin(π * x) * f(y, z) +vy(x, y, z) = sin(π * y) * f(z, x) +vz(x, y, z) = sin(π * z) * f(x, y) +# diagonal deviatoric stress +τxx(x, y, z, η) = 2 * η * π * cos(π * x) * f(y, z) +τyy(x, y, z, η) = 2 * η * π * cos(π * y) * f(z, x) +τzz(x, y, z, η) = 2 * η * π * cos(π * z) * f(x, y) +# off-diagonal deviatoric stress +τxy(x, y, z, η) = 4 * η * cos(π * z) * g(x, y) +τxz(x, y, z, η) = 4 * η * cos(π * y) * g(z, x) +τyz(x, y, z, η) = 4 * η * cos(π * x) * g(y, z) +# forcing terms +ρgx(x, y, z, η) = -2 * η * sin(π * x) * (f(y, z) * π^2 - p(y, z)) +ρgy(x, y, z, η) = -2 * η * sin(π * y) * (f(z, x) * π^2 - p(z, x)) +ρgz(x, y, z, η) = -2 * η * sin(π * z) * (f(x, y) * π^2 - p(x, y)) + +function visme() + xs = LinRange(-1, 1, 201) + ys = LinRange(-1, 1, 201) + zs = LinRange(-1, 1, 201) + + Vxm = [vx(x, y, z) for x in xs, y in ys, z in zs] + Vym = [vy(x, y, z) for x in xs, y in ys, z in zs] + Vzm = [vz(x, y, z) for x in xs, y in ys, z in zs] + + Vmag = sqrt.(Vxm .^ 2 .+ Vym .^ 2 .+ Vzm .^ 2) + Vmag_max = maximum(Vmag) + Vxm .*= 0.2 / Vmag_max + Vym .*= 0.2 / Vmag_max + Vzm .*= 0.2 / Vmag_max + + η0 = 1.0 + + τxxm = [τxx(x, y, z, η0) for x in xs, y in ys, z in zs] + τyym = [τyy(x, y, z, η0) for x in xs, y in ys, z in zs] + τzzm = [τzz(x, y, z, η0) for x in xs, y in ys, z in zs] + τxym = [τxy(x, y, z, η0) for x in xs, y in ys, z in zs] + τxzm = [τxz(x, y, z, η0) for x in xs, y in ys, z in zs] + τyzm = [τyz(x, y, z, η0) for x in xs, y in ys, z in zs] + + τII = sqrt.(0.5 .* (τxxm .^ 2 .+ τyym .^ 2 .+ τzzm .^ 2) .+ τxym .^ 2 .+ τxzm .^ 2 .+ τyzm .^ 2) + + fig = Figure(; size=(500, 450)) + ax = Axis3(fig[1, 1][1, 1]; aspect=:data, xlabel="x", ylabel="y", zlabel="", title="free slip 3D") + limits!(ax, -1, 1, -1, 1, -1, 1) + # hm = heatmap!(ax, xs, ys, τII; colormap=:roma) + vl = volume!(ax, xs, ys, zs, τII; algorithm=:absorption, colormap=:roma) + Colorbar(fig[1, 1][1, 2], vl) + st = 20 + pt = vec([Point3(x, y, z) for x in xs[1:st:end], y in ys[1:st:end], z in zs[1:st:end]]) + dr = Vec3.(vec(Vxm[1:st:end, 1:st:end, 1:st:end]), + vec(Vym[1:st:end, 1:st:end, 1:st:end]), + vec(Vzm[1:st:end, 1:st:end, 1:st:end])) + # arrows!(ax, pt, dr; color=:grey, arrowsize=0.05) + + display(fig) + + # save("free_surface_sol.png", fig) + # save("free_slip_sol.png", fig) + save("free_slip_sol_3D.png", fig) + return +end + +visme() diff --git a/src/Architectures.jl b/src/Architectures.jl deleted file mode 100644 index 9e85210b..00000000 --- a/src/Architectures.jl +++ /dev/null @@ -1,43 +0,0 @@ -module Architectures - -export Architecture - -export launch!, set_device!, get_device, set_device_and_priority!, heuristic_groupsize -export backend, device, details - -using FastIce.Grids - -using KernelAbstractions - -struct Architecture{Kind,B,D,Details} - backend::B - device::D - details::Details -end - -struct SingleDevice end - -function Architecture(backend::Backend, device_id::Integer=1) - device = get_device(backend, device_id) - return Architecture{SingleDevice,typeof(backend),typeof(device),Nothing}(backend, device, nothing) -end - -device(arch::Architecture) = arch.device -backend(arch::Architecture) = arch.backend -details(arch::Architecture) = arch.details - -set_device!(arch::Architecture) = set_device!(arch.device) - -function set_device_and_priority!(arch::Architecture, prio::Symbol) - set_device!(arch) - KernelAbstractions.priority!(arch.backend, prio) - return -end - -set_device!(::Architecture{Kind,CPU}) where {Kind} = nothing -get_device(::CPU, id::Integer) = nothing - -heuristic_groupsize(arch::Architecture, ::Val{N}) where {N} = heuristic_groupsize(arch.device, Val(N)) -heuristic_groupsize(::Architecture{Kind,CPU}, ::Val{N}) where {Kind,N} = 256 - -end diff --git a/src/BoundaryConditions/BoundaryConditions.jl b/src/BoundaryConditions/BoundaryConditions.jl deleted file mode 100644 index c5eaae55..00000000 --- a/src/BoundaryConditions/BoundaryConditions.jl +++ /dev/null @@ -1,63 +0,0 @@ -module BoundaryConditions - -export FieldBoundaryCondition, BoundaryConditionsBatch -export apply_boundary_conditions!, apply_all_boundary_conditions! -export merge_boundary_conditions - -export DirichletBC, HalfCell, FullCell -export ExtrapolateBC, ExtrapolateLogBC, ExpandBC -export ContinuousBC, DiscreteBC -export BoundaryFunction, DiscreteBoundaryFunction, ContinuousBoundaryFunction - -export HideBoundaries, hide - -using FastIce.Grids -using FastIce.Fields -using FastIce.Utils -using FastIce.Architectures - -using KernelAbstractions -using Adapt - -abstract type FieldBoundaryCondition end - -include("field_boundary_conditions.jl") -include("utils.jl") -include("boundary_function.jl") -include("dirichlet_bc.jl") -include("extrapolate_bc.jl") -include("hide_boundaries.jl") - -struct BoundaryConditionsBatch{F,BC} - fields::F - conditions::BC -end - -function merge_boundary_conditions(bc1::BoundaryConditionsBatch, bc2::BoundaryConditionsBatch) - BoundaryConditionsBatch((bc1.fields..., bc2.fields...), - (bc1.conditions..., bc2.conditions...)) -end - -merge_boundary_conditions(bc1::BoundaryConditionsBatch, ::Nothing) = bc1 - -merge_boundary_conditions(::Nothing, bc2::BoundaryConditionsBatch) = bc2 - -@inline function apply_boundary_conditions!(::Val{S}, ::Val{D}, - arch::Architecture, - grid::CartesianGrid, - batch::BoundaryConditionsBatch; kwargs...) where {S,D} - apply_boundary_conditions!(Val(S), Val(D), arch, grid, batch.fields, batch.conditions; kwargs...) -end - -@inline function apply_boundary_conditions!(::Val{S}, ::Val{D}, - arch::Architecture, - grid::CartesianGrid, - batches::NTuple{N,BoundaryConditionsBatch}; kwargs...) where {S,D,N} - ntuple(Val(N)) do I - apply_boundary_conditions!(Val(S), Val(D), arch, grid, batches[I].fields, batches[I].conditions; kwargs...) - end -end - -apply_boundary_conditions!(side, val, arch, grid, ::Nothing; kwargs...) = nothing - -end diff --git a/src/BoundaryConditions/boundary_function.jl b/src/BoundaryConditions/boundary_function.jl deleted file mode 100644 index 0f38ec32..00000000 --- a/src/BoundaryConditions/boundary_function.jl +++ /dev/null @@ -1,39 +0,0 @@ -abstract type BoundaryFunction{F} end - -struct ReducedDimensions end -struct FullDimensions end - -@inline _reduce(::Type{ReducedDimensions}, I, dim) = remove_dim(dim, I) -@inline _reduce(::Type{FullDimensions}, I, dim) = I - -struct ContinuousBoundaryFunction{F,P,RF} <: BoundaryFunction{F} - fun::F - parameters::P - ContinuousBoundaryFunction{RF}(fun::F, params::P) where {RF,F,P} = new{F,P,RF}(fun, params) -end - -struct DiscreteBoundaryFunction{F,P,RF} <: BoundaryFunction{F} - fun::F - parameters::P - DiscreteBoundaryFunction{RF}(fun::F, params::P) where {F,P,RF} = new{F,P,RF}(fun, params) -end - -const CBF{F,P,RF} = ContinuousBoundaryFunction{F,P,RF} where {F,P,RF} -const DBF{F,P,RF} = DiscreteBoundaryFunction{F,P,RF} where {F,P,RF} - -const CDBF{F,P} = Union{CBF{F,P},DBF{F,P}} where {F,P} - -@inline _params(::CDBF{F,Nothing}) where {F} = () -@inline _params(cbf::CDBF{F}) where {F} = cbf.parameters - -@inline (bc::CBF{F,P,RF})(grid, loc, dim, I) where {F,P,RF} = bc.fun(_reduce(RF, coord(grid, loc, I), dim)..., _params(bc)...) - -@inline (bc::DBF{F,P,RF})(grid, loc, dim, I) where {F,P,RF} = bc.fun(grid, loc, dim, Tuple(_reduce(RF, I, dim))..., _params(bc)...) - -# Create a continous or discrete boundary function -# if discrete = true, the function has signature f(grid, loc, dim, inds...) -# if reduce_dims = false, the boundary condition function accepts the same number of coordinates as the number of indices -function BoundaryFunction(fun::Function; discrete=false, parameters=nothing, reduce_dims=true) - RF = reduce_dims ? ReducedDimensions : FullDimensions - discrete ? DiscreteBoundaryFunction{RF}(fun, parameters) : ContinuousBoundaryFunction{RF}(fun, parameters) -end diff --git a/src/BoundaryConditions/dirichlet_bc.jl b/src/BoundaryConditions/dirichlet_bc.jl deleted file mode 100644 index 51b10c7f..00000000 --- a/src/BoundaryConditions/dirichlet_bc.jl +++ /dev/null @@ -1,34 +0,0 @@ -# Reconstruct gradient from the interface between two grid locations -struct HalfCell end - -# Reconstruct gradient from the ghost node location -struct FullCell end - -# First-kind Dirichlet boundary conditon parametrised by the gradient reconstruction kind (can be HalfCell or FullCell) -struct DirichletBC{Gradient,T} <: FieldBoundaryCondition - condition::T - DirichletBC{Gradient}(condition::T) where {Gradient,T} = new{Gradient,T}(condition) -end - -# Create a DirichletBC with a continuous or discrete boundary function -function DirichletBC{G}(fun::Function; kwargs...) where {G} - condition = BoundaryFunction(fun; kwargs...) - return DirichletBC{G}(condition) -end - -@inline (bc::DirichletBC{G,<:Number})(grid, loc, dim, I) where {G} = bc.condition -Base.@propagate_inbounds (bc::DirichletBC{G,<:AbstractArray})(grid, loc, dim, I) where {G} = bc.condition[remove_dim(dim, I)] - -Base.@propagate_inbounds (bc::DirichletBC{G,<:BoundaryFunction})(grid, loc, dim, I) where {G} = bc.condition(grid, loc, dim, I) - -Adapt.adapt_structure(to, f::DirichletBC{G,<:AbstractArray}) where {G} = DirichletBC{G}(Adapt.adapt(to, parent(f.condition))) - -Base.@propagate_inbounds _get_gradient(f2, bc::DirichletBC{HalfCell}, grid, loc, dim, I) = 2 * (bc(grid, loc, dim, I) - f2) -Base.@propagate_inbounds _get_gradient(f2, bc::DirichletBC{FullCell}, grid, loc, dim, I) = (bc(grid, loc, dim, I) - f2) - -@inline function _apply_field_boundary_condition!(side, dim, grid, f, loc, Ibc, bc::DirichletBC) - I = _bc_index(dim, side, loc, size(f), Ibc) - DI = _bc_offset(Val(ndims(grid)), dim, side) - @inbounds f[I] = f[I+DI] + _get_gradient(f[I+DI], bc, grid, loc, dim, I) - return -end diff --git a/src/BoundaryConditions/extrapolate_bc.jl b/src/BoundaryConditions/extrapolate_bc.jl deleted file mode 100644 index a8af4a22..00000000 --- a/src/BoundaryConditions/extrapolate_bc.jl +++ /dev/null @@ -1,29 +0,0 @@ -# Boundary conditions that extrapolates the values of the field outside of the domain -struct ExtrapolateBC <: FieldBoundaryCondition end - -@inline function _apply_field_boundary_condition!(side, dim, grid, f, loc, Ibc, ::ExtrapolateBC) - I = _bc_index(dim, side, loc, size(f), Ibc) - DI = _bc_offset(Val(ndims(grid)), dim, side) - @inbounds f[I] = 2 * f[I+DI] - f[I+2DI] - return -end - -# Boundary conditions that extrapolates the values of the field outside of the domain in logarigthmic space -struct ExtrapolateLogBC <: FieldBoundaryCondition end - -@inline function _apply_field_boundary_condition!(side, dim, grid, f, loc, Ibc, ::ExtrapolateLogBC) - I = _bc_index(dim, side, loc, size(f), Ibc) - DI = _bc_offset(Val(ndims(grid)), dim, side) - @inbounds f[I] = exp(2 * log(f[I+DI]) - log(f[I+2DI])) - return -end - -# Boundary conditions that extrapolates the values of the field outside of the domain by copying them -struct ExpandBC <: FieldBoundaryCondition end - -@inline function _apply_field_boundary_condition!(side, dim, grid, f, loc, Ibc, ::ExpandBC) - I = _bc_index(dim, side, loc, size(f), Ibc) - DI = _bc_offset(Val(ndims(grid)), dim, side) - @inbounds f[I] = f[I+DI] - return -end diff --git a/src/BoundaryConditions/field_boundary_conditions.jl b/src/BoundaryConditions/field_boundary_conditions.jl deleted file mode 100644 index ff287060..00000000 --- a/src/BoundaryConditions/field_boundary_conditions.jl +++ /dev/null @@ -1,44 +0,0 @@ -function apply_boundary_conditions!(::Val{S}, ::Val{D}, - arch::Architecture, - grid::CartesianGrid, - fields::NTuple{N,Field}, - conditions::NTuple{N,FieldBoundaryCondition}; async=true) where {S,D,N} - _validate_fields(fields, D, S) - sizes = ntuple(ifield -> remove_dim(Val(D), size(parent(fields[ifield]))), Val(N)) - halos = ntuple(ifield -> remove_dim(Val(D), halo(fields[ifield])), Val(N)) - offsets = ntuple(Val(N)) do ifield - NF = ndims(fields[ifield]) - 1 - ntuple(Val(NF)) do RD - -halos[ifield][RD][1] - end |> CartesianIndex - end - worksize = max.(sizes...) - _apply_boundary_conditions!(backend(arch), 256)(Val(S), Val(D), grid, sizes, offsets, fields, conditions; ndrange=worksize) - async || synchronize(backend(arch)) - return -end - -@kernel function _apply_boundary_conditions!(side, dim, grid, sizes, offsets, fields, conditions) - I = @index(Global, Cartesian) - # ntuple here unrolls the loop over fields - ntuple(Val(length(fields))) do ifield - Base.@_inline_meta - @inbounds if _checkindices(sizes[ifield], I) - Ibc = I + offsets[ifield] - field = fields[ifield] - condition = conditions[ifield] - _apply_field_boundary_condition!(side, dim, grid, field, location(field), Ibc, condition) - end - end -end - -_apply_field_boundary_condition!(side, dim, grid, field, loc, I, ::Nothing) = nothing - -function _validate_fields(fields::NTuple{N,Field}, dim, side) where {N} - for f in fields - if (location(f, Val(dim)) == Center()) && (halo(f, dim, side) < 1) - error("to apply boundary conditions, halo width must be at least 1") - end - end - return -end diff --git a/src/BoundaryConditions/hide_boundaries.jl b/src/BoundaryConditions/hide_boundaries.jl deleted file mode 100644 index bfc64d79..00000000 --- a/src/BoundaryConditions/hide_boundaries.jl +++ /dev/null @@ -1,30 +0,0 @@ -struct HideBoundaries{N} - workers::NTuple{N,Tuple{Worker,Worker}} - outer_width::NTuple{N,Int} - function HideBoundaries(arch::Architecture, outer_width::NTuple{N,Int}) where {N} - setup() = set_device_and_priority!(arch, :high) - workers = ntuple(D -> (Worker(; setup), Worker(; setup)), Val(N)) - return new{N}(workers, outer_width) - end -end - -function hide(fun::F, hb::HideBoundaries{N}, arch::Architecture, grid::CartesianGrid{N}, boundary_conditions, worksize) where {F,N} - inner_range, outer_ranges = split_ndrange(worksize, hb.outer_width) - # execute inner range in a parent Task with a normal priority - fun(inner_range) - for dim in N:-1:1 - ntuple(Val(2)) do side - worker = hb.workers[dim][side] - range = outer_ranges[dim][side] - batch = boundary_conditions[dim][side] - # execute outer range and boundary conditions asynchronously - put!(worker) do - fun(range) - apply_boundary_conditions!(Val(side), Val(dim), arch, grid, batch) - synchronize(backend(arch)) - end - end - wait.(hb.workers[dim]) # synchronize workers for spatial dimension - end - return -end diff --git a/src/BoundaryConditions/utils.jl b/src/BoundaryConditions/utils.jl deleted file mode 100644 index e498162b..00000000 --- a/src/BoundaryConditions/utils.jl +++ /dev/null @@ -1,23 +0,0 @@ -# Similar to `checkindex`, but for the multidimensional case. Custom axes aren't supported -@inline _checkindices(AI::Tuple, I::Tuple) = (I[1] <= AI[1]) && _checkindices(Base.tail(AI), Base.tail(I)) -@inline _checkindices(AI::Tuple{}, I::Tuple{}) = true -@inline _checkindices(AI::Tuple, I::CartesianIndex) = _checkindices(AI, Tuple(I)) - -# For cell-centered fields, the boundary conditions are specified at index 0 (ghost cell) -@inline _get_i(::Center) = 0 - -# For vertex-centered fields, the boundary conditions are specified at index 1 (first inner point) -@inline _get_i(::Vertex) = 1 - -# Return 1D array index for the "left" and "right" sides of the array -@inline _index1(::Val{1}, L, sz) = _get_i(L) -@inline _index1(::Val{2}, L, sz) = -_get_i(L) + sz + 1 - -# Cartesian array index to store the boundary condition -@inline _bc_index(dim::Val{D}, side, loc, sz, Ibc) where {D} = insert_dim(dim, Ibc, _index1(side, loc[D], sz[D])) - -# Return 1D offset for the "left" and "right" sides of the array to compute the flux projection -@inline _offset1(::Val{1}) = 1 -@inline _offset1(::Val{2}) = -1 - -@inline _bc_offset(N, ::Val{D}, S) where {D} = ntuple(I -> I == D ? _offset1(S) : 0, N) |> CartesianIndex diff --git a/src/Distributed/Distributed.jl b/src/Distributed/Distributed.jl deleted file mode 100644 index a6d42307..00000000 --- a/src/Distributed/Distributed.jl +++ /dev/null @@ -1,38 +0,0 @@ -""" - Distributed - -Tools for performing parallel computations in a distributed environment. -Contains tools for non-blocking halo exchange using MPI. -Implements `BoundaryCondition` API to conveniently define communication as an operation that fills halo buffers. -Together with `HideBoundaries` from the `BoundaryConditions` module enables hiding MPI communication behind computations. -""" -module Distributed - -export CartesianTopology -export global_rank, shared_rank, node_name, cartesian_communicator, shared_communicator, coordinates -export dimensions, global_size, node_size -export global_grid_size, local_grid -export neighbors, neighbor, has_neighbor -export gather! - -export ExchangeInfo, DistributedBoundaryConditions -export DistributedMPI - -using FastIce.Grids -using FastIce.Fields -using FastIce.Architectures -using FastIce.BoundaryConditions -import FastIce.BoundaryConditions: apply_boundary_conditions! - -using MPI -using KernelAbstractions -import ImplicitGlobalGrid - -"Trait structure used as a type parameter to indicate that the Architecture is a distributed MPI Architecture." -struct DistributedMPI end - -include("topology.jl") -include("boundary_conditions.jl") -include("gather.jl") - -end diff --git a/src/Distributed/boundary_conditions.jl b/src/Distributed/boundary_conditions.jl deleted file mode 100644 index d6ad1445..00000000 --- a/src/Distributed/boundary_conditions.jl +++ /dev/null @@ -1,112 +0,0 @@ -""" - ExchangeInfo - -Structure containing the information to exchange halos of one side of an N-dimensional array. -""" -mutable struct ExchangeInfo{SB,RB} - send_buffer::SB - recv_buffer::RB - send_request::MPI.Request - recv_request::MPI.Request - ExchangeInfo(send_buf, recv_buf) = new{typeof(send_buf),typeof(recv_buf)}(send_buf, recv_buf, MPI.REQUEST_NULL, MPI.REQUEST_NULL) -end - -""" - ExchangeInfo(::Val{S}, ::Val{D}, field::Field) where {S,D} - -Create `ExchangeInfo` to exchange halos on side `S` in the spatial direction `D`. -""" -function ExchangeInfo(::Val{S}, ::Val{D}, field::Field) where {S,D} - send_view = get_send_view(Val(S), Val(D), field) - recv_view = get_recv_view(Val(S), Val(D), field) - send_buffer = similar(parent(send_view), eltype(send_view), size(send_view)) - recv_buffer = similar(parent(recv_view), eltype(recv_view), size(recv_view)) - return ExchangeInfo(send_buffer, recv_buffer) -end - -""" - apply_boundary_conditions!(::Val{S}, ::Val{D}, - arch::Architecture, - grid::CartesianGrid, - fields::NTuple{N,Field}, - exchange_infos::NTuple{N,ExchangeInfo}; async=true) where {S,D,N} - -Perform a non-blocking MPI exchange for a set of fields. -""" -function apply_boundary_conditions!(::Val{S}, ::Val{D}, - arch::Architecture, - grid::CartesianGrid, - fields::NTuple{N,Field}, - exchange_infos::NTuple{N,ExchangeInfo}; async=true) where {S,D,N} - comm = cartesian_communicator(details(arch)) - nbrank = neighbor(details(arch), D, S) - - # initiate non-blocking MPI recieve and device-to-device copy to the send buffer - for idx in eachindex(fields) - info = exchange_infos[idx] - info.recv_request = MPI.Irecv!(info.recv_buffer, comm; source=nbrank) - send_view = get_send_view(Val(S), Val(D), fields[idx]) - copyto!(info.send_buffer, send_view) - end - KernelAbstractions.synchronize(backend(arch)) - - # initiate non-blocking MPI send - for idx in eachindex(fields) - info = exchange_infos[idx] - info.send_request = MPI.Isend(info.send_buffer, comm; dest=nbrank) - end - - recv_ready = falses(N) - send_ready = falses(N) - - # test send and receive requests, initiating device-to-device copy - # to the receive buffer if the receive is complete - while !(all(recv_ready) && all(send_ready)) - for idx in eachindex(fields) - info = exchange_infos[idx] - if MPI.Test(info.recv_request) && !recv_ready[idx] - recv_view = get_recv_view(Val(S), Val(D), fields[idx]) - copyto!(recv_view, info.recv_buffer) - recv_ready[idx] = true - end - send_ready[idx] = MPI.Test(info.send_request) - end - yield() - end - async || KernelAbstractions.synchronize(backend(arch)) - - return -end - -_overlap(::Vertex) = 1 -_overlap(::Center) = 0 - -get_recv_view(side::Val{S}, dim::Val{D}, f::Field) where {S,D} = get_recv_view(side, dim, parent(f), halo(f, D, S)) - -function get_send_view(side::Val{S}, dim::Val{D}, f::Field) where {S,D} - get_send_view(side, dim, parent(f), halo(f, D, S), _overlap(location(f, dim))) -end - -function get_recv_view(::Val{1}, ::Val{D}, array::AbstractArray, halo_width::Integer) where {D} - recv_range = Base.OneTo(halo_width) - indices = ntuple(I -> I == D ? recv_range : Colon(), Val(ndims(array))) - return view(array, indices...) -end - -function get_recv_view(::Val{2}, ::Val{D}, array::AbstractArray, halo_width::Integer) where {D} - recv_range = (size(array, D)-halo_width+1):size(array, D) - indices = ntuple(I -> I == D ? recv_range : Colon(), Val(ndims(array))) - return view(array, indices...) -end - -function get_send_view(::Val{1}, ::Val{D}, array::AbstractArray, halo_width::Integer, overlap::Integer) where {D} - send_range = (overlap+halo_width+1):(overlap+2halo_width) - indices = ntuple(I -> I == D ? send_range : Colon(), Val(ndims(array))) - return view(array, indices...) -end - -function get_send_view(::Val{2}, ::Val{D}, array::AbstractArray, halo_width::Integer, overlap::Integer) where {D} - send_range = (size(array, D)-overlap-2halo_width+1):(size(array, D)-overlap-halo_width) - indices = ntuple(I -> I == D ? send_range : Colon(), Val(ndims(array))) - return view(array, indices...) -end diff --git a/src/Distributed/gather.jl b/src/Distributed/gather.jl deleted file mode 100644 index b9b8150e..00000000 --- a/src/Distributed/gather.jl +++ /dev/null @@ -1,20 +0,0 @@ -""" - gather!(dst::Union{AbstractArray{T,N},Nothing}, src::AbstractArray{T,N}, comm::MPI.Comm; root=0) where {T,N} - -Gather local array `src` into a global array `dst`. -Size of the global array `size(dst)` should be equal to the product of the size of a local array `size(src)` and the dimensions of a Cartesian communicator `comm`. -The array will be gathered on the process with id `root` (`root=0` by default). -Note that the memory for a global array should be allocated only on the process with id `root`, on other processes `dst` can be set to `nothing`. -""" -function gather!(dst::Union{AbstractArray{T,N},Nothing}, src::AbstractArray{T,N}, comm::MPI.Comm; root=0) where {T,N} - ImplicitGlobalGrid.gather!(src, dst, comm; root=root) -end - -""" - gather!(arch::Architecture{DistributedMPI}, dst, src::Field; kwargs...) - -Gather the interior of a field `src` into a global array `dst`. -""" -function gather!(arch::Architecture{DistributedMPI}, dst, src::Field; kwargs...) - gather!(dst, interior(src), cartesian_communicator(details(arch)); kwargs...) -end diff --git a/src/Distributed/topology.jl b/src/Distributed/topology.jl deleted file mode 100644 index 025ac5fe..00000000 --- a/src/Distributed/topology.jl +++ /dev/null @@ -1,157 +0,0 @@ -""" - CartesianTopology - -Represents N-dimensional Cartesian topology of distributed processes. -""" -struct CartesianTopology{N} - nprocs::Int - dims::NTuple{N,Int} - global_rank::Int - shared_rank::Int - cart_coords::NTuple{N,Int} - neighbors::NTuple{N,NTuple{2,Int}} - comm::MPI.Comm - cart_comm::MPI.Comm - shared_comm::MPI.Comm - node_name::String -end - -""" - CartesianTopology(dims::NTuple{N,Int}; comm=MPI.COMM_WORLD) where {N} - -Create an N-dimensional Cartesian topology with dimensions `dims` and using base MPI communicator `comm` (by default, `comm=MPI.COMM_WORLD`). -If all entries in `dims` are not equal to `0`, the product of `dims` should be equal to the total number of MPI processes `MPI.Comm_size(comm)`. -If any (or all) entries of `dims` are `0`, the dimensions in the corresponding spatial directions will be picked automatically. -""" -function CartesianTopology(dims::NTuple{N,Int}; comm=MPI.COMM_WORLD) where {N} - nprocs = MPI.Comm_size(comm) - dims = Tuple(MPI.Dims_create(nprocs, dims)) - cart_comm = MPI.Cart_create(MPI.COMM_WORLD, dims) - global_rank = MPI.Comm_rank(cart_comm) - shared_comm = MPI.Comm_split_type(cart_comm, MPI.COMM_TYPE_SHARED, global_rank) - shared_rank = MPI.Comm_rank(shared_comm) - node_name = MPI.Get_processor_name() - cart_coords = Tuple(MPI.Cart_coords(cart_comm)) - - neighbors = ntuple(Val(N)) do dim - MPI.Cart_shift(cart_comm, dim - 1, 1) - end - - return CartesianTopology{N}(nprocs, dims, global_rank, shared_rank, cart_coords, neighbors, comm, cart_comm, shared_comm, node_name) -end - -""" - global_rank(t::CartesianTopology) - -Global id of a process in a Cartesian topology. -""" -global_rank(t::CartesianTopology) = t.global_rank - -""" - shared_rank(t::CartesianTopology) - -Local id of a process within a single node. Can be used to set the GPU device. -""" -shared_rank(t::CartesianTopology) = t.shared_rank - -""" - node_name(t::CartesianTopology) - -Name of a node according to `MPI.Get_processor_name()`. -""" -node_name(t::CartesianTopology) = t.node_name - -""" - cartesian_communicator(t::CartesianTopology) - -MPI Cartesian communicator for the topology. -""" -cartesian_communicator(t::CartesianTopology) = t.cart_comm - -""" - shared_communicator(t::CartesianTopology) - -MPI communicator for the processes sharing the same node. -""" -shared_communicator(t::CartesianTopology) = t.shared_comm - -""" - dimensions(t::CartesianTopology) - -Dimensions of the topology as NTuple. -""" -dimensions(t::CartesianTopology) = t.dims - -""" - coordinates(t::CartesianTopology) - -Coordinates of a current process within a Cartesian topology. -""" -coordinates(t::CartesianTopology) = t.cart_coords - -""" - neighbors(t::CartesianTopology) - -Neighbors of a current process. - -Returns NTuple containing process ids of the two immediate neighbors in each spatial direction, or MPI.PROC_NULL if no neighbor on a corresponding side. -""" -neighbors(t::CartesianTopology) = t.neighbors - -""" - neighbor(t::CartesianTopology, dim, side) - -Returns id of a neighbor process in spatial direction `dim` on the side `side`, if this neighbor exists, or MPI.PROC_NULL otherwise. -""" -neighbor(t::CartesianTopology, dim, side) = t.neighbors[dim][side] - -""" - has_neighbor(t::CartesianTopology, dim, side) - -Returns true if there a neighbor process in spatial direction `dim` on the side `side`, or false otherwise. -""" -has_neighbor(t::CartesianTopology, dim, side) = t.neighbors[dim][side] != MPI.PROC_NULL - -""" - global_size(t::CartesianTopology) - -Total number of processes withing the topology. -""" -global_size(t::CartesianTopology) = MPI.Comm_size(t.cart_comm) - -""" - node_size(t::CartesianTopology) - -Number of processes sharing the same node. -""" -node_size(t::CartesianTopology) = MPI.Comm_size(t.shared_comm) - -""" - global_grid_size(t::CartesianTopology, local_size) - -Return the global size for a structured grid. -""" -global_grid_size(t::CartesianTopology, local_size) = t.dims .* local_size - -""" - local_grid(g::CartesianGrid, t::CartesianTopology) - -Return a `CartesianGrid` covering the subdomain which corresponds to the current process. -""" -function local_grid(g::CartesianGrid, t::CartesianTopology) - local_extent = Tuple(extent(g)) ./ t.dims - local_size = size(g) .÷ t.dims - local_origin = Tuple(origin(g)) .+ local_extent .* t.cart_coords - - return CartesianGrid(local_origin, local_extent, local_size) -end - -""" - Architecture(backend::Backend, topology::CartesianTopology) where {N} - -Create a distributed Architecture using `backend` and `topology`. For GPU backends, device will be selected automatically based on a process id within a node. -""" -function Architectures.Architecture(backend::Backend, topology::CartesianTopology) - device = get_device(backend, shared_rank(topology) + 1) - return Architecture{DistributedMPI,typeof(backend),typeof(device),typeof(topology)}(backend, device, topology) -end diff --git a/src/FastIce.jl b/src/FastIce.jl index 40de0aa9..06149306 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -25,19 +25,12 @@ https://github.com/PTsolvers/FastIce.jl greet(; kwargs...) = printstyled(GREETING; kwargs...) greet_fast(; kwargs...) = printstyled(GREETING_FAST; kwargs...) -# core modules -include("Grids/Grids.jl") -include("GridOperators.jl") +# core modules (included in alphabetical order) +include("Geometries/Geometries.jl") +include("LevelSets/LevelSets.jl") include("Logging.jl") -include("Architectures.jl") -include("Fields/Fields.jl") -include("Utils/Utils.jl") -include("BoundaryConditions/BoundaryConditions.jl") -include("KernelLaunch.jl") -include("Distributed/Distributed.jl") include("Physics.jl") include("Writers.jl") -include("LevelSets/LevelSets.jl") # ice flow models include("Models/Models.jl") diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl deleted file mode 100644 index ccb0adbb..00000000 --- a/src/Fields/Fields.jl +++ /dev/null @@ -1,156 +0,0 @@ -module Fields - -export AbstractField -export Field, interior -export FunctionField -export location, data, halo, set! -export ConstantField, ZeroField, OneField - -using FastIce.Grids -using FastIce.GridOperators -using FastIce.Architectures - -using Adapt -using OffsetArrays -using KernelAbstractions - -import LinearAlgebra - -abstract type AbstractField{T,N,L} <: AbstractArray{T,N} end - -Base.@pure location(::AbstractField{T,N,L}) where {T,N,L} = L.instance -Base.@pure location(::AbstractField{T,N,L}, ::Val{I}) where {T,N,L,I} = L.instance[I] -Base.@pure Base.eltype(::AbstractField{T}) where {T} = T -Base.@pure Base.ndims(::AbstractField{T,N}) where {T,N} = N - -Base.IndexStyle(::AbstractField) = IndexCartesian() - -struct Field{T,N,L,D,H,I} <: AbstractField{T,N,L} - data::D - halo::H - indices::I - Field{L}(data::D, halo::H, indices::I) where {L,D,H,I} = new{eltype(data),ndims(data),L,D,H,I}(data, halo, indices) -end - -data_axis(sz, h) = (1-h[1]):(sz+h[2]) - -function make_data(backend, T, sz, halo_sz) - total_halo_size = map(sum, halo_sz) - array = KernelAbstractions.allocate(backend, T, sz .+ total_halo_size) - field_axes = ntuple(I -> data_axis(sz[I], halo_sz[I]), Val(length(sz))) - return OffsetArray(array, field_axes) -end - -const HaloSize{N,I<:Integer} = NTuple{N,Tuple{I,I}} - -function Field(backend::Backend, grid::CartesianGrid, T::DataType, loc::L, halo::HaloSize) where {L} - sz = size(grid, loc) - data = make_data(backend, T, sz, halo) - indices = Base.OneTo.(sz) - return Field{L}(data, halo, indices) -end - -Field(arch::Architecture, args...; kwargs...) = Field(backend(arch), args...; kwargs...) - -expand_axis_halo(::Nothing) = (0, 0) -expand_axis_halo(halo::Integer) = (halo, halo) -expand_axis_halo(halo::Tuple) = (something(halo[1], 0), something(halo[2], 0)) - -expand_halo(::Val{N}, halo::HaloSize{N}) where {N} = halo -expand_halo(::Val{N}, halo::Tuple) where {N} = ntuple(I -> expand_axis_halo(halo[I]), Val(length(halo))) -expand_halo(::Val{N}, halo::Integer) where {N} = ntuple(I -> (halo, halo), Val(N)) -expand_halo(::Val{N}, halo::Nothing) where {N} = ntuple(I -> (0, 0), Val(N)) - -expand_loc(::Val{N}, loc::NTuple{N,Location}) where {N} = loc -expand_loc(::Val{N}, loc::Location) where {N} = ntuple(_ -> loc, Val(N)) - -function Field(backend::Backend, grid::CartesianGrid, loc::L, T::DataType=eltype(grid); halo::H=nothing) where {L,H} - N = ndims(grid) - return Field(backend, grid, T, expand_loc(Val(N), loc), expand_halo(Val(N), halo)) -end - -Base.checkbounds(f::Field, I...) = checkbounds(f.data, I...) -Base.checkbounds(f::Field, I::Union{CartesianIndex,AbstractArray{<:CartesianIndex}}) = checkbounds(f.data, I) - -Base.checkbounds(::Type{Bool}, f::Field, I...) = checkbounds(Bool, f.data, I...) -Base.checkbounds(::Type{Bool}, f::Field, I::Union{CartesianIndex,AbstractArray{<:CartesianIndex}}) = checkbounds(Bool, f.data, I) - -Base.size(f::Field) = length.(f.indices) -Base.parent(f::Field) = parent(f.data) -Base.axes(f::Field) = f.indices - -Base.view(f::Field, I...) = view(f.data, I...) - -data(f::Field) = f.data - -halo(f::Field) = f.halo -halo(f::Field, dim::Integer) = f.halo[dim] -halo(f::Field, dim::Integer, side::Integer) = f.halo[dim][side] - -Adapt.adapt_structure(to, f::Field{T,N,L}) where {T,N,L} = Field{L}(Adapt.adapt(to, f.data), f.halo, f.indices) - -Base.@propagate_inbounds Base.getindex(f::Field, inds...) = getindex(f.data, inds...) -Base.@propagate_inbounds Base.setindex!(f::Field, val, inds...) = setindex!(f.data, val, inds...) - -Base.@propagate_inbounds Base.firstindex(f::Field) = firstindex(f.data) -Base.@propagate_inbounds Base.firstindex(f::Field, dim) = firstindex(f.data, dim) -Base.@propagate_inbounds Base.lastindex(f::Field) = lastindex(f.data) -Base.@propagate_inbounds Base.lastindex(f::Field, dim) = lastindex(f.data, dim) - -function interior_indices(f::Field) - return ntuple(ndims(f)) do I - (firstindex(parent(f), I)+f.halo[I][1]):(lastindex(parent(f), I)-f.halo[I][2]) - end -end - -function interior(f::Field) - return view(parent(f), interior_indices(f)...) -end - -set!(f::Field, other::Field) = (copy!(interior(f), interior(other)); nothing) -set!(f::Field{T}, val::T) where {T<:Number} = (fill!(interior(f), val); nothing) -set!(f::Field, A::AbstractArray) = (copy!(interior(f), A); nothing) - -@kernel function _set_continuous!(dst, grid, loc, fun::F, args...) where {F} - I = @index(Global, Cartesian) - dst[I] = fun(coord(grid, loc, I)..., args...) -end - -@kernel function _set_discrete!(dst, grid, loc, fun::F, args...) where {F} - I = @index(Global, Cartesian) - dst[I] = fun(grid, loc, I, args...) -end - -function set!(f::Field{T,N}, grid::CartesianGrid{N}, fun::F; discrete=false, parameters=()) where {T,F,N} - loc = location(f) - dst = interior(f) - if discrete - _set_discrete!(get_backend(dst), 256, size(dst))(dst, grid, loc, fun, parameters...) - else - _set_continuous!(get_backend(dst), 256, size(dst))(dst, grid, loc, fun, parameters...) - end - return -end - -# grid operators -Base.@propagate_inbounds ∂(f::AbstractField, I, dim) = ∂(f, I, dim, location(f, dim)) -Base.@propagate_inbounds ∂(f::AbstractField, I, dim, ::Vertex) = ∂ᶜ(f, I, dim) -Base.@propagate_inbounds ∂(f::AbstractField, I, dim, ::Center) = ∂ᵛ(f, I, dim) - -# field norm -LinearAlgebra.norm(f::Field) = LinearAlgebra.norm(interior(f)) -LinearAlgebra.norm(f::Field, p::Real) = LinearAlgebra.norm(interior(f), p) - -function Base.show(io::IO, ::MIME"text/plain", field::Field{T,N,L}) where {T,N,L} - print(io, "$(N)D $(join(size(field), "×")) Field{$T}\n") -end - -function Base.show(io::IO, field::Field{T,N,L}) where {T,N,L} - print(io, "$(N)D $(join(size(field), "×")) Field{$T}") -end - -include("function_field.jl") - -include("constant_field.jl") - -end diff --git a/src/Fields/constant_field.jl b/src/Fields/constant_field.jl deleted file mode 100644 index 45d257ba..00000000 --- a/src/Fields/constant_field.jl +++ /dev/null @@ -1,16 +0,0 @@ -struct ZeroField{T} <: AbstractField{T,1,Nothing} end - -@inline Base.getindex(::ZeroField{T}, inds...) where {T} = zero(T) -@inline Base.size(::ZeroField) = (1,) - -struct OneField{T} <: AbstractField{T,1,Nothing} end - -@inline Base.getindex(::OneField{T}, inds...) where {T} = one(T) -@inline Base.size(::OneField) = (1,) - -struct ConstantField{T} <: AbstractField{T,1,Nothing} - value::T -end - -@inline Base.getindex(f::ConstantField, inds...) = f.value -@inline Base.size(::ConstantField) = (1,) diff --git a/src/Fields/function_field.jl b/src/Fields/function_field.jl deleted file mode 100644 index e522583e..00000000 --- a/src/Fields/function_field.jl +++ /dev/null @@ -1,46 +0,0 @@ -struct Continuous end -struct Discrete end - -struct FunctionField{T,N,L,CD,F,G,P} <: AbstractField{T,N,L} - func::F - grid::G - parameters::P - function FunctionField{CD,L}(func::F, grid::G, parameters::P) where {CD,L,F,G,P} - N = ndims(grid) - T = eltype(grid) - return new{T,N,L,CD,F,G,P}(func, grid, parameters) - end -end - -function FunctionField(func::F, grid::G, loc; discrete=false, parameters=nothing) where {F,G} - loc = expand_loc(Val(ndims(grid)), loc) - L = typeof(loc) - CD = discrete ? Discrete : Continuous - return FunctionField{CD,L}(func, grid, parameters) -end - -Base.size(f::FunctionField) = size(f.grid, location(f)) - -@inline func_type(::FunctionField{T,N,L,CD}) where {T,N,L,CD} = CD - -@inline _params(::Nothing) = () -@inline _params(p) = p - -Base.@propagate_inbounds function call_func(func::F, grid, loc, I, params, ::Type{Continuous}) where {F} - func(coord(grid, loc, I)..., params...) -end - -Base.@propagate_inbounds function call_func(func::F, grid, loc, I, params, ::Type{Discrete}) where {F} - func(grid, loc, I, params...) -end - -Base.@propagate_inbounds Base.getindex(f::FunctionField{T,N}, inds::Vararg{Integer,N}) where {T,N} = getindex(f, CartesianIndex(inds)) -Base.@propagate_inbounds function Base.getindex(f::FunctionField{T,N}, I::CartesianIndex{N}) where {T,N} - call_func(f.func, f.grid, location(f), I, _params(f.parameters), func_type(f)) -end - -function Adapt.adapt_structure(to, f::FunctionField{T,N,L,CD}) where {T,N,L,CD} - FunctionField{CD,L}(Adapt.adapt(to, f.func), - Adapt.adapt(to, f.grid), - Adapt.adapt(to, f.parameters)) -end diff --git a/src/Geometries/Geometries.jl b/src/Geometries/Geometries.jl new file mode 100644 index 00000000..9b964afb --- /dev/null +++ b/src/Geometries/Geometries.jl @@ -0,0 +1,15 @@ +module Geometries + +# export AbstractElevation, AABB, DataElevation, SyntheticElevation +# export domain, rotated_domain, rotation, generate_elevation, extents, center, dilate +export SyntheticElevation +export extents, dilate +export make_synthetic + +using Chmy.Grids + +include("elevation_data.jl") +include("geometry_helpers.jl") +include("synthetic_geometries.jl") + +end # module diff --git a/src/Geometries/elevation_data.jl b/src/Geometries/elevation_data.jl new file mode 100644 index 00000000..847bb050 --- /dev/null +++ b/src/Geometries/elevation_data.jl @@ -0,0 +1,62 @@ +struct AABB{T<:Union{Real}} + xmin::T + xmax::T + ymin::T + ymax::T + zmin::T + zmax::T +end + +""" + AABB(xs, ys, zs) + +Construct an axis aligend bounding box `AABB` from the coordinates `xs`, `ys`, and `zs`. +""" +AABB(xs, ys, zs) = AABB(extrema(xs)..., extrema(ys)..., extrema(zs)...) + +abstract type AbstractElevation{T<:Real} end + +struct DataElevation{T,V<:AbstractArray{T},R<:AbstractRange{T},M<:AbstractMatrix{T}} <: AbstractElevation{T} + x::R + y::R + offsets::V + z_bed::M + z_surf::M + rotation::M + domain::AABB{T} + rotated_domain::AABB{T} +end + +"Generate DataElevation data." +function DataElevation(x, y, offsets, z_bed, z_surf, R) + # get non-rotated domain + domain = AABB(extrema(x)..., extrema(y)..., Float64(minimum([minimum(filtered(z_bed)), minimum(filtered(z_surf))])), + Float64(maximum([maximum(filtered(z_bed)), maximum(filtered(z_surf))]))) + # rotate bed and surface + bed_extents = AABB(rotate_minmax(x, y, z_bed, R)...) + surf_extents = AABB(rotate_minmax(x, y, z_surf, R)...) + # get rotated domain + rotated_domain = union(bed_extents, surf_extents) + dem = DataElevation(x, y, offsets, z_bed, z_surf, R, domain, rotated_domain) + return dem +end + +struct SyntheticElevation{T,B,S} <: AbstractElevation{T} + z_bed::B + z_surf::S + domain::AABB{T} +end + +"Generate SyntheticElevation data." +function SyntheticElevation(lx, ly, zmin, zmax, z_bed, z_surf) + domain = AABB(-lx / 2, lx / 2, -ly / 2, ly / 2, zmin, zmax) + return SyntheticElevation(z_bed, z_surf, domain) +end + +domain(dem::DataElevation) = dem.domain +rotated_domain(dem::DataElevation) = dem.rotated_domain +rotation(dem::DataElevation) = dem.rotation + +domain(dem::SyntheticElevation) = dem.domain +rotated_domain(dem::AbstractElevation) = domain(dem) +rotation(::AbstractElevation) = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] diff --git a/src/Geometries/geometry_helpers.jl b/src/Geometries/geometry_helpers.jl new file mode 100644 index 00000000..2c4e7033 --- /dev/null +++ b/src/Geometries/geometry_helpers.jl @@ -0,0 +1,39 @@ +""" + extents(box::AABB) + +Returns the extents of the bounding box `box`. +""" +extents(box::AABB) = box.xmax - box.xmin, box.ymax - box.ymin, box.zmax - box.zmin + +"AABB center" +function center(box::AABB{T}) where {T} + half = convert(T, 0.5) + return half * (box.xmin + box.xmax), half * (box.ymin + box.ymax), half * (box.zmin + box.zmax) +end + +"Dilate AABB by extending its limits around the center by certain fraction in each dimension" +function dilate(box::AABB, fractions) + Δx, Δy, Δz = extents(box) .* fractions + return AABB(box.xmin - Δx, box.xmax + Δx, box.ymin - Δy, box.ymax + Δy, box.zmin - Δz, box.zmax + Δz) +end + +"Filter NaNs." +filtered(X) = filter(_x -> !isnan(_x), X) + +"Create AABB enclosing both box1 and box2" +function union(box1::AABB, box2::AABB) + return AABB(min(box1.xmin, box2.xmin), max(box1.xmax, box2.xmax), + min(box1.ymin, box2.ymin), max(box1.ymax, box2.ymax), + min(box1.zmin, box2.zmin), max(box1.zmax, box2.zmax)) +end + +"Rotate field `X`, `Y`, `Z` with rotation matrix `R`." +function rotate(X, Y, Z, R) + xrot = R[1, 1] .* X .+ R[1, 2] .* Y .+ R[1, 3] .* Z + yrot = R[2, 1] .* X .+ R[2, 2] .* Y .+ R[2, 3] .* Z + zrot = R[3, 1] .* X .+ R[3, 2] .* Y .+ R[3, 3] .* Z + return xrot, yrot, zrot +end + +"Rotate field `X`, `Y`, `Z` with rotation matrix `R` and return extents." +rotate_minmax(X, Y, Z, R) = rotate(collect(extrema(X)), collect(extrema(Y)), collect(extrema(filtered(Z))), R) diff --git a/src/Geometries/synthetic_geometries.jl b/src/Geometries/synthetic_geometries.jl new file mode 100644 index 00000000..58c41956 --- /dev/null +++ b/src/Geometries/synthetic_geometries.jl @@ -0,0 +1,89 @@ +function generate_turtle(nx, ny, lx, ly, zmin, zmax) + @info "expects lx, ly, lz = 5.0, 5.0, 1.0" + # Numerics + amp = 0.1 + ω = 10π + tanβ = tan(-π/6) + el = 0.35 + gl = 0.9 + x = LinRange(-lx / 2, lx / 2, nx + 1) + y = LinRange(-lx / 2, ly / 2, ny + 1) + # Functions + generate_bed(x, y) = amp * sin(ω * x) * sin(ω * y) + tanβ * x + el + (1.5 * y)^2 + generate_surf(x, y) = 1.0 - ((1.7 * x + 0.22)^2 + (0.5 * y)^2) + # Generate + bed = [(zmax - zmin) * generate_bed(x / lx, y / ly) for x in x, y in y] + surf = [gl * generate_surf(x / lx, y / ly) for x in x, y in y] + return bed, surf +end + +function generate_dome(nx, ny, lx, ly) + # Numerics + h_0 = 0.8 + x = LinRange(-lx / 2, lx / 2, nx + 1) + y = LinRange(-lx / 2, ly / 2, ny + 1) + # Functions + generate_bed(x, y) = 0.0 + generate_surf(x, y) = h_0 - (x^2 + y^2) + # Compute + bed = [generate_bed(_x, _y) for _x in x, _y in y] + surf = [generate_surf(_x, _y) for _x in x, _y in y] + surf[surf. i == D ? op(I[i], 1) : I[i], Val(N)) |> CartesianIndex -end - -Base.@assume_effects :foldable function δ(op, I::CartesianIndex{N}, ::Val{D1}, ::Val{D2}) where {N,D1,D2} - δI = ntuple(Val(N)) do i - (i == D1 || i == D2) ? op(I[i], 1) : I[i] - end - return CartesianIndex(δI) -end - -@propagate_inbounds ∂ᶜ(fv, I, D) = fv[δ(+, I, D)] - fv[I] -@propagate_inbounds ∂ᵛ(fc, I, D) = fc[I] - fc[δ(-, I, D)] - -@propagate_inbounds avᶜ(fv, I, D) = 0.5 * (fv[δ(+, I, D)] + fv[I]) -@propagate_inbounds avᵛ(fc, I, D) = 0.5 * (fc[I] + fc[δ(-, I, D)]) - -@propagate_inbounds avᶜ(fv, I, D1, D2) = 0.25 * (fv[I] + fv[δ(+, I, D1)] + fv[δ(+, I, D2)] + fv[δ(+, I, D1, D2)]) -@propagate_inbounds avᵛ(fc, I, D1, D2) = 0.25 * (fc[I] + fc[δ(-, I, D1)] + fc[δ(-, I, D2)] + fc[δ(-, I, D1, D2)]) - -@propagate_inbounds maxlᶜ(fv, I, D) = max(fv[δ(+, I, D)], fv[I]) -@propagate_inbounds maxlᵛ(fc, I, D) = max(fc[I], fc[δ(-, I, D)]) - -@propagate_inbounds maxlᵛ(fc, I, D1, D2) = max(fc[δ(-, I, D2)], fc[I], fc[δ(+, I, D2)], - fc[δ(-, I, D1, D2)], fc[δ(-, I, D1)], fc[δ(+, δ(-, I, D1), D2)]) - -for (dim, val) in ((:x, 1), (:y, 2), (:z, 3)) - for loc in (:ᶜ, :ᵛ) - ∂l = Symbol(:∂, loc) - avl = Symbol(:av, loc) - maxll = Symbol(:maxl, loc) - - ∂ = Symbol(∂l, dim) - av = Symbol(avl, dim) - maxl = Symbol(maxll, dim) - - @eval begin - @propagate_inbounds $∂(f, I) = $∂l(f, I, Val($val)) - @propagate_inbounds $av(f, I) = $avl(f, I, Val($val)) - @propagate_inbounds $maxl(f, I) = $maxll(f, I, Val($val)) - end - end -end - -for (dim, val1, val2) in ((:xy, 1, 2), (:xz, 1, 3), (:yz, 2, 3)) - for loc in (:ᶜ, :ᵛ) - avl = Symbol(:av, loc) - av = Symbol(avl, dim) - - @eval begin - @propagate_inbounds $av(f, I) = $avl(f, I, Val($val1), Val($val2)) - end - end -end - -for (dim, val1, val2) in ((:xy, 1, 2), (:xz, 1, 3), (:yz, 2, 3), (:yx, 2, 1), (:zx, 3, 1), (:zy, 3, 2)) - for loc in (:ᶜ, :ᵛ) - maxll = Symbol(:maxl, loc) - maxl = Symbol(maxll, dim) - - @eval @propagate_inbounds $maxl(f, I) = $maxll(f, I, Val($val1), Val($val2)) - end -end - -end diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl deleted file mode 100644 index e37785fc..00000000 --- a/src/Grids/Grids.jl +++ /dev/null @@ -1,18 +0,0 @@ -module Grids - -export Location, Center, Vertex -export DiscreteAxis, spacing, Δ, origin, extent, coord, center, vertex, coords, centers, vertices - -export CartesianGrid, axis -export xcoord, ycoord, zcoord, xcenter, ycenter, zcenter, xvertex, yvertex, zvertex -export xcoords, ycoords, zcoords, xcenters, ycenters, zcenters, xvertices, yvertices, zvertices - -abstract type Location end - -struct Center <: Location end -struct Vertex <: Location end - -include("discrete_axis.jl") -include("cartesian_grid.jl") - -end diff --git a/src/Grids/cartesian_grid.jl b/src/Grids/cartesian_grid.jl deleted file mode 100644 index 9a86cd6c..00000000 --- a/src/Grids/cartesian_grid.jl +++ /dev/null @@ -1,193 +0,0 @@ -""" -Rectilinear grid with uniform spacing. - -Examples -======== - -```jldoctest -julia> using FastIce.Grids - -julia> grid = CartesianGrid(origin=(0.0,0.0), extent=(1.0,1.0), size=(4,4)) -2D 4×4 CartesianGrid{Float64}: - x ∈ [0.0–1.0]; Δx = 0.25 - y ∈ [0.0–1.0]; Δy = 0.25 -``` -""" -struct CartesianGrid{N,T<:AbstractFloat,I<:Integer} - axes::NTuple{N,DiscreteAxis{T,I}} -end - -function CartesianGrid(origin::NTuple{N,T}, extent::NTuple{N,T}, size::NTuple{N,I}) where {N,T,I} - CartesianGrid(DiscreteAxis.(origin, extent, size)) -end - -""" - CartesianGrid(; origin::NTuple{N,T}, extent::NTuple{N,T}, size::NTuple{N,I}) - -Create a Cartesian grid with a specified origin (bottom-south-west corner in 3D), spatial extent, and number of grid cells. -""" -CartesianGrid(; origin::NTuple{N,T}, extent::NTuple{N,T}, size::NTuple{N,I}) where {N,T,I} = CartesianGrid(origin, extent, size) - -# overload methods from Base -@pure Base.eltype(::CartesianGrid{N,T}) where {N,T} = T -@pure Base.ndims(::CartesianGrid{N}) where {N} = N - -Base.size(grid::CartesianGrid) = length.(grid.axes) - -@propagate_inbounds Base.size(grid::CartesianGrid, dim::Integer) = length(grid.axes[dim]) -@propagate_inbounds Base.size(grid::CartesianGrid, loc::Location, dim::Integer) = length(grid.axes[dim], loc) - -Base.size(grid::CartesianGrid{N}, locs::NTuple{N,Location}) where {N} = length.(grid.axes, locs) -Base.size(grid::CartesianGrid, loc::Location) = ntuple(D -> length(grid.axes[D], loc), Val(ndims(grid))) - -# pretty-printing -function Base.show(io::IO, grid::CartesianGrid{N,T}) where {N,T} - print(io, "$(N)D $(join(size(grid), "×")) CartesianGrid{$T}:\n") - symbols = (:x, :y, :z) - for idim in 1:N - l, r = origin(grid, idim), origin(grid, idim) + extent(grid, idim) - print(io, " $(symbols[idim]) ∈ [$(l)–$(r)]; Δ$(symbols[idim]) = $(spacing(grid, idim))\n") - end -end - -_name_coords(tp::NTuple{1}) = NamedTuple{(:x, )}(tp) -_name_coords(tp::NTuple{2}) = NamedTuple{(:x, :y)}(tp) -_name_coords(tp::NTuple{3}) = NamedTuple{(:x, :y, :z)}(tp) - -""" - axis(grid::CartesianGrid, dim::Integer) - -Return a `DiscreteAxis` corresponding to the spatial dimension `dim`. -""" -axis(grid::CartesianGrid, dim::Integer) = grid.axes[dim] - -""" - origin(grid::CartesianGrid) - -Return the origin of a `CartesianGrid`. The origin corresponds to bottom-south-west corner of the grid in 3D. -""" -origin(grid::CartesianGrid) = _name_coords(origin.(grid.axes)) - -""" - extent(grid::CartesianGrid) - -Return the spatial extent of a `CartesianGrid` as a tuple. -""" -extent(grid::CartesianGrid) = _name_coords(extent.(grid.axes)) - -""" - spacing(grid::CartesianGrid) - -Return the spacing of a `CartesianGrid` as a tuple. -""" -spacing(grid::CartesianGrid) = _name_coords(spacing.(grid.axes)) - -""" - Δ(grid::CartesianGrid) - -Return the spacing of a `CartesianGrid` as a tuple. Same as `spacing`. -""" -Δ(grid::CartesianGrid) = spacing(grid) - -@propagate_inbounds origin(grid::CartesianGrid, dim::Integer) = origin(grid.axes[dim]) -@propagate_inbounds extent(grid::CartesianGrid, dim::Integer) = extent(grid.axes[dim]) -@propagate_inbounds spacing(grid::CartesianGrid, dim::Integer) = spacing(grid.axes[dim]) -@propagate_inbounds Δ(grid::CartesianGrid, dim::Integer) = spacing(grid.axes[dim]) - -""" - coord(grid::CartesianGrid{N}, loc::NTuple{N,Location}, inds::NTuple{N}) - -Return a tuple of spatial coordinates of a grid point at location `loc` and indices `inds`. - -For vertex-centered locations, first grid point is at the origin. -For cell-centered locations, first grid point at half-spacing distance from the origin. -""" -@propagate_inbounds function coord(grid::CartesianGrid{N}, loc::NTuple{N,Location}, inds::NTuple{N}) where {N} - ntuple(Val(N)) do I - Base.@_inline_meta - coord(grid.axes[I], loc[I], inds[I]) - end -end - -@propagate_inbounds coord(grid::CartesianGrid{N}, loc::Location, inds::NTuple{N}) where {N} = coord.(grid.axes, Ref(loc), inds) - -coord(grid::CartesianGrid{N}, loc, I::CartesianIndex{N}) where {N} = coord(grid, loc, Tuple(I)) - -@propagate_inbounds coord(grid::CartesianGrid, loc::Location, dim::Integer, i::Integer) = coord(grid.axes[dim], loc, i) -@propagate_inbounds function coord(grid::CartesianGrid{N}, loc::Location, dim::Integer, inds::NTuple{N}) where {N} - coord(grid.axes[dim], loc, inds[dim]) -end -@propagate_inbounds function coord(grid::CartesianGrid{N}, locs::NTuple{N,Location}, dim::Integer, inds::NTuple{N}) where {N} - coord(grid.axes[dim], locs[dim], inds[dim]) -end -@propagate_inbounds function coord(grid::CartesianGrid{N}, locs::NTuple{N,Location}, dim::Integer, i::Integer) where {N} - coord(grid.axes[dim], locs[dim], i) -end - -@propagate_inbounds coord(grid::CartesianGrid, loc, ::Val{D}, i) where {D} = coord(grid, loc, D, i) - -center(grid::CartesianGrid{N}, inds::NTuple{N}) where {N} = center.(grid.axes, inds) -vertex(grid::CartesianGrid{N}, inds::NTuple{N}) where {N} = vertex.(grid.axes, inds) - -center(grid::CartesianGrid{N}, I::CartesianIndex{N}) where {N} = center(grid, Tuple(I)) -vertex(grid::CartesianGrid{N}, I::CartesianIndex{N}) where {N} = vertex(grid, Tuple(I)) - -@propagate_inbounds center(grid::CartesianGrid, dim::Integer, i) = coord(grid, Center(), dim, i) -@propagate_inbounds center(grid::CartesianGrid, ::Val{D}, i) where {D} = coord(grid, Center(), D, i) - -@propagate_inbounds vertex(grid::CartesianGrid, dim::Integer, i) = coord(grid, Vertex(), dim, i) -@propagate_inbounds vertex(grid::CartesianGrid, ::Val{D}, i) where {D} = coord(grid, Vertex(), D, i) - -xcoord(grid::CartesianGrid, loc, i) = coord(grid, loc, Val(1), i) -ycoord(grid::CartesianGrid, loc, i) = coord(grid, loc, Val(2), i) -zcoord(grid::CartesianGrid, loc, i) = coord(grid, loc, Val(3), i) - -xcenter(grid::CartesianGrid, i) = center(grid, Val(1), i) -ycenter(grid::CartesianGrid, i) = center(grid, Val(2), i) -zcenter(grid::CartesianGrid, i) = center(grid, Val(3), i) - -xvertex(grid::CartesianGrid, i) = vertex(grid, Val(1), i) -yvertex(grid::CartesianGrid, i) = vertex(grid, Val(2), i) -zvertex(grid::CartesianGrid, i) = vertex(grid, Val(3), i) - -coords(grid::CartesianGrid, loc::Location, dim::Integer; kwargs...) = coords(grid.axes[dim], loc; kwargs...) -coords(grid::CartesianGrid, loc::Location, ::Val{D}; kwargs...) where {D} = coords(grid.axes[D], loc; kwargs...) - -@inline function coords(grid::CartesianGrid{N}, locs::NTuple{N,Location}; halos=nothing) where {N} - ntuple(Val(N)) do I - Base.@_inline_meta - isnothing(halos) ? coords(grid, locs[I], Val(I)) : coords(grid, locs[I], Val(I); halo=halos[I]) - end -end - -xcoords(grid::CartesianGrid, loc::Location; kwargs...) = coords(grid, loc, Val(1); kwargs...) -ycoords(grid::CartesianGrid, loc::Location; kwargs...) = coords(grid, loc, Val(2); kwargs...) -zcoords(grid::CartesianGrid, loc::Location; kwargs...) = coords(grid, loc, Val(3); kwargs...) - -centers(grid::CartesianGrid, dim::Integer; kwargs...) = centers(grid.axes[dim]; kwargs...) -centers(grid::CartesianGrid, ::Val{D}; kwargs...) where {D} = centers(grid.axes[D]; kwargs...) - -@inline function centers(grid::CartesianGrid{N}; halos=nothing) where {N} - ntuple(Val(N)) do I - Base.@_inline_meta - isnothing(halos) ? centers(grid, Val(I)) : centers(grid, Val(I); halo=halos[I]) - end -end - -xcenters(grid::CartesianGrid; kwargs...) = centers(grid, Val(1); kwargs...) -ycenters(grid::CartesianGrid; kwargs...) = centers(grid, Val(2); kwargs...) -zcenters(grid::CartesianGrid; kwargs...) = centers(grid, Val(3); kwargs...) - -vertices(grid::CartesianGrid, dim::Integer; kwargs...) = vertices(grid.axes[dim]; kwargs...) -vertices(grid::CartesianGrid, ::Val{D}; kwargs...) where {D} = vertices(grid.axes[D]; kwargs...) - -@inline function vertices(grid::CartesianGrid{N}; halos=nothing) where {N} - ntuple(Val(N)) do I - Base.@_inline_meta - isnothing(halos) ? vertices(grid, Val(I)) : vertices(grid, Val(I); halo=halos[I]) - end -end - -xvertices(grid::CartesianGrid; kwargs...) = vertices(grid, Val(1); kwargs...) -yvertices(grid::CartesianGrid; kwargs...) = vertices(grid, Val(2); kwargs...) -zvertices(grid::CartesianGrid; kwargs...) = vertices(grid, Val(3); kwargs...) diff --git a/src/Grids/discrete_axis.jl b/src/Grids/discrete_axis.jl deleted file mode 100644 index 87e1e728..00000000 --- a/src/Grids/discrete_axis.jl +++ /dev/null @@ -1,61 +0,0 @@ -""" -Discretised one-dimensional interval. Grid type `CartesianGrid` is defined as a tuple of `DiscreteAxis`'s. -""" -struct DiscreteAxis{T<:AbstractFloat,I<:Integer} <: AbstractRange{T} - origin::T - extent::T - step::T - len::I - DiscreteAxis(origin::T, extent::T, len::I) where {T,I} = new{T,I}(origin, extent, extent / len, len) -end - -# overload methods from Base - -import Base.@pure -import Base.@propagate_inbounds - -@pure Base.eltype(::DiscreteAxis{T}) where {T} = T - -Base.length(ax::DiscreteAxis) = ax.len -Base.length(ax::DiscreteAxis, ::Center) = ax.len -Base.length(ax::DiscreteAxis{T,I}, ::Vertex) where {T,I} = ax.len + oneunit(I) - -Base.step(ax::DiscreteAxis) = ax.step - -@propagate_inbounds Base.getindex(ax::DiscreteAxis{T}, i::Integer) where {T} = ax.origin + ax.step / 2 + (i - 1) * ax.step - -spacing(ax::DiscreteAxis) = ax.step -Δ(ax::DiscreteAxis) = ax.step - -origin(ax::DiscreteAxis) = ax.origin -origin(ax::DiscreteAxis, ::Vertex) = ax.origin -origin(ax::DiscreteAxis, ::Center) = ax.origin + ax.step / 2 - -extent(ax::DiscreteAxis) = ax.extent -extent(ax::DiscreteAxis, ::Vertex) = ax.extent -extent(ax::DiscreteAxis, ::Center) = ax.extent - ax.step - -coord(ax::DiscreteAxis{T}, loc::Location, i::Integer) where {T} = origin(ax, loc) + (i - 1) * ax.step -center(ax, i::Integer) = coord(ax, Center(), i) -vertex(ax, i::Integer) = coord(ax, Vertex(), i) - -@inline function coords(ax::DiscreteAxis, loc::Location; halo=nothing) - if isnothing(halo) - halo = (0, 0) - end - - if halo isa Integer - halo = (halo, halo) - end - - return _coords(ax, loc, halo) -end - -@inline function _coords(ax::DiscreteAxis, loc::Location, halo::Tuple{Integer,Integer}) - start = coord(ax, loc, 1 - halo[1]) - stop = coord(ax, loc, length(ax, loc) + halo[2]) - LinRange(start, stop, length(ax, loc) + sum(halo)) -end - -centers(ax::DiscreteAxis; kwargs...) = coords(ax, Center(); kwargs...) -vertices(ax::DiscreteAxis; kwargs...) = coords(ax, Vertex(); kwargs...) diff --git a/src/KernelLaunch.jl b/src/KernelLaunch.jl deleted file mode 100644 index 498f2cd5..00000000 --- a/src/KernelLaunch.jl +++ /dev/null @@ -1,86 +0,0 @@ -module KernelLaunch - -export launch! - -using FastIce.Grids -using FastIce.Architectures -using FastIce.BoundaryConditions - -using KernelAbstractions - -""" - launch!(arch::Architecture, grid::CartesianGrid, kernel::Pair{K,Args}; ) where {K,Args} - -Launch a KernelAbstraction `kernel` on a `grid` using the backend from `arch`. -Either `worksize` or `location` must be provided as keyword arguments. - -# Keyword Arguments -- `worksize`: worksize of a kernel, i.e. how many grid points are included in each spatial direction. -- `location[=nothing]`: compute worksize as a size of the grid at a specified location. - If only one location is provided, e.g. `location=Vertex()`, then this location will be used for all spacial directions. -- `offset[=nothing]`: index offset for all grid indices as a `CartesianIndex`. -- `expand[=nothing]`: if provided, the worksize is increased by `2*expand`, and offset is set to `-expand`, or combined with user-provided offset. -- `hide_boundaries[=nothing]`: instance of `HideBoundaries`, that will be used to overlap boundary processing with computations at inner points of the domain. -- `outer_width[=nothing]`: if `hide_boundaries` is specified, used to determine the decomposition of the domain into inner and outer regions. -- `boundary_conditions[=nothing]`: a tuple of boundary condition batches for each side of every spatial direction. -- `async[=false]`: if set to `true`, will return immediately, without waiting on the kernel execution to finish. -""" -function launch!(arch::Architecture, grid::CartesianGrid, kernel::Pair{K,Args}; - worksize=nothing, - location=nothing, - offset=nothing, - expand=nothing, - boundary_conditions=nothing, - hide_boundaries=nothing, - async=false) where {K,Args} - fun, args = kernel - - if isnothing(location) && isnothing(worksize) - throw(ArgumentError("either grid location or worksize must me specified")) - end - - if !isnothing(location) && !isnothing(worksize) - throw(ArgumentError("either grid location or worksize must me specified, but not both")) - end - - if isnothing(worksize) - worksize = size(grid, location) - end - - if !isnothing(expand) - if !isnothing(offset) - offset -= oneunit(CartesianIndex{ndims(grid)}) - else - offset = -oneunit(CartesianIndex{ndims(grid)}) - end - worksize = worksize .+ 2 .* expand - end - - groupsize = heuristic_groupsize(arch, Val(length(worksize))) - - if isnothing(hide_boundaries) - fun(backend(arch), groupsize, worksize)(args..., offset) - isnothing(boundary_conditions) || apply_all_boundary_conditions!(arch, grid, boundary_conditions) - else - hide(hide_boundaries, arch, grid, boundary_conditions, worksize) do indices - sub_offset, ndrange = first(indices) - oneunit(first(indices)), size(indices) - if !isnothing(offset) - sub_offset += offset - end - fun(backend(arch), groupsize)(args..., sub_offset; ndrange) - end - end - - async || KernelAbstractions.synchronize(backend(arch)) - return -end - -function apply_all_boundary_conditions!(arch::Architecture, grid::CartesianGrid{N}, boundary_conditions) where {N} - ntuple(Val(N)) do D - apply_boundary_conditions!(Val(1), Val(D), arch, grid, boundary_conditions[D][1]) - apply_boundary_conditions!(Val(2), Val(D), arch, grid, boundary_conditions[D][2]) - end - return -end - -end diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index d399aced..fca590e2 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -1,15 +1,22 @@ module LevelSets -export compute_level_set_from_dem! +export compute_levelset_from_dem!, invert_levelset! +export ω_from_ψ! + +using Chmy.Architectures +using Chmy.BoundaryConditions +using Chmy.Distributed +using Chmy.Grids +using Chmy.Fields +using Chmy.GridOperators +using Chmy.KernelLaunch -using FastIce -using FastIce.Grids -using FastIce.Architectures -using FastIce.Fields using KernelAbstractions using LinearAlgebra, GeometryBasics include("signed_distances.jl") -include("compute_level_sets.jl") +include("compute_levelsets.jl") +include("volume_fractions.jl") +include("volume_fractions_kernels.jl") -end \ No newline at end of file +end diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl deleted file mode 100644 index ccd2d26b..00000000 --- a/src/LevelSets/compute_level_sets.jl +++ /dev/null @@ -1,24 +0,0 @@ -""" - _init_level_set!(Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, cutoff::AbstractFloat, R::AbstractMatrix) - -Initialize level sets. -""" -@kernel function _init_level_set!(Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, cutoff, R) - I = @index(Global, Cartesian) - x, y, z = coord(Ψ_grid, location(Ψ), I) - P = R * Point3(x, y, z) - ud, sgn = sd_dem(P, cutoff, dem, dem_grid) - Ψ[I] = ud * sgn -end - -""" - compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, R=LinearAlgebra.I) - -Compute level sets from dem. -""" -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, R=LinearAlgebra.I) - cutoff = 4maximum(spacing(Ψ_grid)) - kernel = _init_level_set!(backend(arch), 256, size(Ψ)) - kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) - return -end diff --git a/src/LevelSets/compute_levelsets.jl b/src/LevelSets/compute_levelsets.jl new file mode 100644 index 00000000..8a41ab0f --- /dev/null +++ b/src/LevelSets/compute_levelsets.jl @@ -0,0 +1,35 @@ +@kernel inbounds = true function init_levelset!(Ψ::Field, dem::AbstractField, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R, O) + I = @index(Global, NTuple) + I = I + O + x, y, z = coord(Ψ_grid, location(Ψ), I...) + P = R * Point3(x, y, z) + ud, sgn = sd_dem(P, cutoff, dem, dem_grid) + Ψ[I...] = ud * sgn +end + +@kernel inbounds = true function invert_levelset!(Ψ::Field, O) + I = @index(Global, NTuple) + I = I + O + Ψ[I...] = -Ψ[I...] +end + +""" + compute_levelset_from_dem!(arch::Architecture, launch, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) + +Compute level sets from dem. +""" +function compute_levelset_from_dem!(arch::Architecture, launch, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid{2}, grid::UniformGrid{3}, R=LinearAlgebra.I) + cutoff = 4maximum(spacing(grid)) + launch(arch, grid, init_levelset! => (Ψ, dem, dem_grid2D, grid, cutoff, R); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) + return +end + +""" + invert_levelset!(arch::Architecture, launch, Ψ::Field, grid::UniformGrid) + +Invert level set `Ψ` to set what's below the surface as "inside". +""" +function invert_levelset!(arch::Architecture, launch, Ψ::Field, grid::UniformGrid) + launch(arch, grid, invert_levelset! => (Ψ,); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) + return +end diff --git a/src/LevelSets/signed_distances.jl b/src/LevelSets/signed_distances.jl index 5e4b9da6..8ff43e85 100644 --- a/src/LevelSets/signed_distances.jl +++ b/src/LevelSets/signed_distances.jl @@ -1,7 +1,26 @@ -export sd_dem +# 2D functions +function sd_poly(p::Point2{T}, poly::AbstractVector{Point2{T}}) where {T} + d = dot(p - poly[1], p - poly[1]) + s = 1.0 + j = length(poly) + for i in eachindex(poly) + e = poly[j] - poly[i] + w = p - poly[i] + b = w - e .* clamp(dot(w, e) / dot(e, e), zero(T), oneunit(T)) + d = min(d, dot(b, b)) + c = p[2] >= poly[i][2], p[2] < poly[j][2], e[1] * w[2] > e[2] * w[1] + if all(c) || all(.!c) + s = -s + end + j = i + end + return s * sqrt(d) +end + +# 3D functions @inline S(x) = x == zero(x) ? oneunit(x) : sign(x) -@inline sign_triangle(p, a, b, c) = S(dot(p - a, cross(b - a, c - a))) +sign_triangle(p, a, b, c) = S(dot(p - a, cross(b - a, c - a))) @inline function ud_triangle(p, a, b, c) dot2(v) = dot(v, v) @@ -12,23 +31,22 @@ export sd_dem ac = a - c pc = p - c nor = cross(ba, ac) - return sqrt( - (sign(dot(cross(ba, nor), pa)) + - sign(dot(cross(cb, nor), pb)) + - sign(dot(cross(ac, nor), pc)) < 2) - ? - min( - dot2(ba * clamp(dot(ba, pa) / dot2(ba), 0, 1) - pa), - dot2(cb * clamp(dot(cb, pb) / dot2(cb), 0, 1) - pb), - dot2(ac * clamp(dot(ac, pc) / dot2(ac), 0, 1) - pc)) - : - dot(nor, pa) * dot(nor, pa) / dot2(nor)) + return sqrt((sign(dot(cross(ba, nor), pa)) + + sign(dot(cross(cb, nor), pb)) + + sign(dot(cross(ac, nor), pc)) < 2) + ? + min(dot2(ba * clamp(dot(ba, pa) / dot2(ba), 0, 1) - pa), + dot2(cb * clamp(dot(cb, pb) / dot2(cb), 0, 1) - pb), + dot2(ac * clamp(dot(ac, pc) / dot2(ac), 0, 1) - pc)) + : + dot(nor, pa) * dot(nor, pa) / dot2(nor)) end @inline function closest_vertex_index(P, grid) - Δ = spacing(grid) - O = Grids.origin(grid) - I = @. Int(fld(P - O, Δ)) + 1 + Δs = spacing(grid, Vertex(), 1, 1) + O = Grids.origin(grid, Vertex()) + X = @. fld(P - O, Δs) + I = @. Int(X) + 1 I1 = 1 I2 = size(grid, Vertex()) return clamp.(I, I1, I2) |> CartesianIndex @@ -53,9 +71,11 @@ end return ud, sign_triangle(P, T_BL...) end -Base.clamp(p::NTuple{N}, grid::CartesianGrid{N}) where {N} = clamp.(p, Grids.origin(grid), Grids.origin(grid) .+ extent(grid)) +function Base.clamp(p::NTuple{N}, grid::UniformGrid{N}) where {N} + clamp.(p, Grids.origin(grid, Vertex()), Grids.origin(grid, Vertex()) .+ extent(grid, Vertex())) +end -function sd_dem(P, cutoff, dem, grid) +function sd_dem(P, cutoff, dem, grid::UniformGrid{2}) @inbounds Pp = clamp((P[1], P[2]), grid) BL = closest_vertex_index(Pp .- cutoff, grid) TR = closest_vertex_index(Pp .+ cutoff, grid) @@ -69,4 +89,4 @@ function sd_dem(P, cutoff, dem, grid) ud = min(ud, ud_pair) end return ud, sgn -end \ No newline at end of file +end diff --git a/src/LevelSets/volume_fractions.jl b/src/LevelSets/volume_fractions.jl new file mode 100644 index 00000000..0ebc7077 --- /dev/null +++ b/src/LevelSets/volume_fractions.jl @@ -0,0 +1,127 @@ +@inline perturb(ψ) = abs(ψ) > oftype(ψ, 1e-20) ? ψ : (ψ > zero(ψ) ? oftype(ψ, 1e-20) : oftype(ψ, -1e-20)) + +# volume of 2- and 3-simplex +@inline vol(v::Vararg{Vec2{T},3}) where {T} = T(inv(2)) * abs(det([v[3] - v[1] v[2] - v[1]])) +@inline vol(v::Vararg{Vec3{T},4}) where {T} = T(inv(6)) * abs(det([v[2] - v[1] v[3] - v[1] v[4] - v[1]])) + +# volume of 2-simplex with immersed boundary given by signed distances +function vol(ψ::NTuple{3,T}, v::NTuple{3,Vec2{T}}) where {T} + if ψ[1] < 0 && ψ[2] < 0 && ψ[3] < 0 # --- + return vol(v[1], v[2], v[3]) + elseif ψ[1] > 0 && ψ[2] > 0 && ψ[3] > 0 # +++ + return zero(T) + end + @inline vij(i, j) = v[j] * (ψ[i] / (ψ[i] - ψ[j])) - v[i] * (ψ[j] / (ψ[i] - ψ[j])) + v12, v13, v23 = vij(1, 2), vij(1, 3), vij(2, 3) + if ψ[1] < 0 + if ψ[2] < 0 + vol(v[1], v23, v13) + vol(v[1], v[2], v23) # --+ + else + if ψ[3] < 0 + vol(v[3], v12, v23) + vol(v[3], v[1], v12) # -+- + else + vol(v[1], v12, v13) # -++ + end + end + else + if ψ[2] < 0 + if ψ[3] < 0 + vol(v[2], v13, v12) + vol(v[2], v[3], v13) # +-- + else + vol(v12, v[2], v23) # +-+ + end + else + vol(v13, v23, v[3]) # ++- + end + end +end + +# volume of 2d box with immersed boundary given by signed distances +function vol(ψ::NTuple{4,T}, Δ::NTuple{2,T}) where {T} + v00, v01, v10, v11 = Vec2(zero(T), zero(T)), Vec2(Δ[1], zero(T)), Vec2(zero(T), Δ[2]), Vec2(Δ[1], Δ[2]) + ψ00, ψ01, ψ10, ψ11 = perturb.(ψ) + return vol((ψ00, ψ01, ψ11), (v00, v01, v11)) + + vol((ψ00, ψ11, ψ10), (v00, v11, v10)) +end + +# volume of 3-simplex with immersed boundary given by signed distances +function vol(ψ::NTuple{4,T}, v::NTuple{4,Vec3{T}}) where {T} + @inline vij(i, j) = v[j] * (ψ[i] / (ψ[i] - ψ[j])) - v[i] * (ψ[j] / (ψ[i] - ψ[j])) + nneg = count(ψ .< 0) + if nneg == 0 # ++++ + return zero(T) + elseif nneg == 1 # -+++ + if ψ[1] < 0 + return vol(v[1], vij(1, 2), vij(1, 3), vij(1, 4)) + elseif ψ[2] < 0 + return vol(v[2], vij(2, 1), vij(2, 3), vij(2, 4)) + elseif ψ[3] < 0 + return vol(v[3], vij(3, 1), vij(3, 2), vij(3, 4)) + else # ψ[4] < 0 + return vol(v[4], vij(4, 1), vij(4, 2), vij(4, 3)) + end + elseif nneg == 2 # --++ + if ψ[1] < 0 && ψ[2] < 0 + return vol(v[1], v[2], vij(1, 3), vij(2, 4)) + + vol(vij(2, 3), v[2], vij(1, 3), vij(2, 4)) + + vol(v[1], vij(1, 4), vij(1, 3), vij(2, 4)) + elseif ψ[1] < 0 && ψ[3] < 0 + return vol(v[1], v[3], vij(1, 4), vij(3, 2)) + + vol(vij(3, 4), v[3], vij(1, 4), vij(3, 2)) + + vol(v[1], vij(1, 2), vij(1, 4), vij(3, 2)) + elseif ψ[1] < 0 && ψ[4] < 0 + return vol(v[1], v[4], vij(1, 2), vij(4, 3)) + + vol(vij(4, 2), v[4], vij(1, 2), vij(4, 3)) + + vol(v[1], vij(1, 3), vij(1, 2), vij(4, 3)) + elseif ψ[2] < 0 && ψ[3] < 0 + return vol(v[3], v[2], vij(3, 1), vij(2, 4)) + + vol(vij(2, 1), v[2], vij(3, 1), vij(2, 4)) + + vol(v[3], vij(3, 4), vij(3, 1), vij(2, 4)) + elseif ψ[2] < 0 && ψ[4] < 0 + return vol(v[4], v[2], vij(4, 1), vij(2, 3)) + + vol(vij(2, 1), v[2], vij(4, 1), vij(2, 3)) + + vol(v[4], vij(4, 3), vij(4, 1), vij(2, 3)) + else # ψ[3] < 0 && ψ[4] < 0 + return vol(v[3], v[4], vij(3, 1), vij(4, 2)) + + vol(vij(4, 1), v[4], vij(3, 1), vij(4, 2)) + + vol(v[3], vij(3, 2), vij(3, 1), vij(4, 2)) + end + elseif nneg == 3 # ---+ + vol_tot = vol(v[1], v[2], v[3], v[4]) + if ψ[1] >= 0 + return vol_tot - vol(v[1], vij(1, 2), vij(1, 3), vij(1, 4)) + elseif ψ[2] >= 0 + return vol_tot - vol(v[2], vij(2, 1), vij(2, 3), vij(2, 4)) + elseif ψ[3] >= 0 + return vol_tot - vol(v[3], vij(3, 1), vij(3, 2), vij(3, 4)) + else # ψ[4] >= 0 + return vol_tot - vol(v[4], vij(4, 1), vij(4, 2), vij(4, 3)) + end + else # ---- + return vol(v[1], v[2], v[3], v[4]) + end +end + +# volume of 3d-box with immersed boundary given by signed distances +function vol(ψ::NTuple{8,T}, Δ::NTuple{3,T}) where {T} + v000, v001, v010, v011 = Vec3(zero(T), zero(T), zero(T)), Vec3(Δ[1], zero(T), zero(T)), Vec3(zero(T), Δ[2], zero(T)), Vec3(Δ[1], Δ[2], zero(T)) + v100, v101, v110, v111 = Vec3(zero(T), zero(T), Δ[3]), Vec3(Δ[1], zero(T), Δ[3]), Vec3(zero(T), Δ[2], Δ[3]), Vec3(Δ[1], Δ[2], Δ[3]) + ψ = perturb.(ψ) + return vol((ψ[1], ψ[5], ψ[3], ψ[2]), (v000, v100, v010, v001)) + + vol((ψ[7], ψ[5], ψ[3], ψ[8]), (v110, v100, v010, v111)) + + vol((ψ[6], ψ[5], ψ[8], ψ[2]), (v101, v100, v111, v001)) + + vol((ψ[4], ψ[8], ψ[3], ψ[2]), (v011, v111, v010, v001)) + + vol((ψ[8], ψ[5], ψ[3], ψ[2]), (v111, v100, v010, v001)) +end + +# volume fraction inside a 2D grid cell +function ψ2ω(ψ::NTuple{2,NTuple{2,T}}, Δ::NTuple{2,T}) where {T} + ψf = (ψ[1][1], ψ[1][2], ψ[2][1], ψ[2][2]) + return vol(ψf, Δ) / prod(Δ) +end + +# volume fraction inside a 3D grid cell +function ψ2ω(ψ::NTuple{2,NTuple{2,NTuple{2,T}}}, Δ::NTuple{3,T}) where {T} + ψf = (ψ[1][1][1], ψ[1][1][2], ψ[1][2][1], ψ[1][2][2], ψ[2][1][1], ψ[2][1][2], ψ[2][2][1], ψ[2][2][2]) + return vol(ψf, Δ) / prod(Δ) +end diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl new file mode 100644 index 00000000..89a723e1 --- /dev/null +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -0,0 +1,110 @@ +# 2D case +@kernel inbounds = true function _ω_from_ψ!(ω, ψ, g::StructuredGrid{2}, O=Offset()) + I = @index(Global, NTuple) + i, j = I + O + + # cell center + cell_cc = ((ψ[i, j], ψ[i+1, j]), (ψ[i, j+1], ψ[i+1, j+1])) + ω.cc[i, j] = ψ2ω(cell_cc, Δ(g, Center(), i, j)) + + # vertex + cell_vv = ((lerp(ψ, Center(), g, i - 1, j - 1), lerp(ψ, Center(), g, i, j - 1)), + (lerp(ψ, Center(), g, i - 1, j + 0), lerp(ψ, Center(), g, i, j + 0))) + ω.vv[i, j] = ψ2ω(cell_vv, Δ(g, Vertex(), i, j)) + + # locations + loc_vc = (Vertex(), Center()) + loc_cv = (Center(), Vertex()) + + # x interface + cell_vc = ((lerp(ψ, loc_cv, g, i - 1, j + 0), lerp(ψ, loc_cv, g, i, j + 0)), + (lerp(ψ, loc_cv, g, i - 1, j + 1), lerp(ψ, loc_cv, g, i, j + 1))) + ω.vc[i, j] = ψ2ω(cell_vc, Δ(g, loc_vc, i, j)) + + # y interface + cell_cv = ((lerp(ψ, loc_vc, g, i, j - 1), lerp(ψ, loc_vc, g, i + 1, j - 1)), + (lerp(ψ, loc_vc, g, i, j + 0), lerp(ψ, loc_vc, g, i + 1, j + 0))) + ω.cv[i, j] = ψ2ω(cell_cv, Δ(g, loc_cv, i, j)) +end + +# 3D case +@kernel inbounds = true function _ω_from_ψ!(ω, ψ, g::StructuredGrid{3}, O=Offset()) + I = @index(Global, NTuple) + i, j, k = I + O + + # cell center + cell_ccc = (((ψ[i, j, k+0], ψ[i+1, j, k+0]), (ψ[i, j+1, k+0], ψ[i+1, j+1, k+0])), + ((ψ[i, j, k+1], ψ[i+1, j, k+1]), (ψ[i, j+1, k+1], ψ[i+1, j+1, k+1]))) + ω.ccc[i, j, k] = ψ2ω(cell_ccc, Δ(g, Center(), i, j, k)) + + # vertex + cell_vvv = (((lerp(ψ, Center(), g, i - 1, j - 1, k - 1), lerp(ψ, Center(), g, i, j - 1, k - 1)), + (lerp(ψ, Center(), g, i - 1, j + 0, k - 1), lerp(ψ, Center(), g, i, j + 0, k - 1))), + ((lerp(ψ, Center(), g, i - 1, j - 1, k + 0), lerp(ψ, Center(), g, i, j - 1, k + 0)), + (lerp(ψ, Center(), g, i - 1, j + 0, k + 0), lerp(ψ, Center(), g, i, j + 0, k + 0)))) + ω.vvv[i, j, k] = ψ2ω(cell_vvv, Δ(g, Vertex(), i, j, k)) + + # locations + loc_vcc = (Vertex(), Center(), Center()) + loc_cvc = (Center(), Vertex(), Center()) + loc_ccv = (Center(), Center(), Vertex()) + + loc_vvc = (Vertex(), Vertex(), Center()) + loc_vcv = (Vertex(), Center(), Vertex()) + loc_cvv = (Center(), Vertex(), Vertex()) + + # x interface + cell_vcc = (((lerp(ψ, loc_cvv, g, i - 1, j + 0, k + 0), lerp(ψ, loc_cvv, g, i, j + 0, k + 0)), + (lerp(ψ, loc_cvv, g, i - 1, j + 1, k + 0), lerp(ψ, loc_cvv, g, i, j + 1, k + 0))), + ((lerp(ψ, loc_cvv, g, i - 1, j + 0, k + 1), lerp(ψ, loc_cvv, g, i, j + 0, k + 1)), + (lerp(ψ, loc_cvv, g, i - 1, j + 1, k + 1), lerp(ψ, loc_cvv, g, i, j + 1, k + 1)))) + ω.vcc[i, j, k] = ψ2ω(cell_vcc, Δ(g, loc_vcc, i, j, k)) + + # y interface + cell_cvc = (((lerp(ψ, loc_vcv, g, i, j - 1, k + 0), lerp(ψ, loc_vcv, g, i + 1, j - 1, k + 0)), + (lerp(ψ, loc_vcv, g, i, j + 0, k + 0), lerp(ψ, loc_vcv, g, i + 1, j + 0, k + 0))), + ((lerp(ψ, loc_vcv, g, i, j - 1, k + 1), lerp(ψ, loc_vcv, g, i + 1, j - 1, k + 1)), + (lerp(ψ, loc_vcv, g, i, j + 0, k + 1), lerp(ψ, loc_vcv, g, i + 1, j + 0, k + 1)))) + ω.cvc[i, j, k] = ψ2ω(cell_cvc, Δ(g, loc_cvc, i, j, k)) + + # z interface + cell_ccv = (((lerp(ψ, loc_vvc, g, i, j + 0, k - 1), lerp(ψ, loc_vvc, g, i + 1, j + 0, k - 1)), + (lerp(ψ, loc_vvc, g, i, j + 1, k - 1), lerp(ψ, loc_vvc, g, i + 1, j + 1, k - 1))), + ((lerp(ψ, loc_vvc, g, i, j + 0, k + 0), lerp(ψ, loc_vvc, g, i + 1, j + 0, k + 0)), + (lerp(ψ, loc_vvc, g, i, j + 1, k + 0), lerp(ψ, loc_vvc, g, i + 1, j + 1, k + 0)))) + ω.ccv[i, j, k] = ψ2ω(cell_ccv, Δ(g, loc_ccv, i, j, k)) + + # xy edge + cell_vvc = (((lerp(ψ, loc_ccv, g, i - 1, j - 1, k + 0), lerp(ψ, loc_ccv, g, i, j - 1, k + 0)), + (lerp(ψ, loc_ccv, g, i - 1, j + 0, k + 0), lerp(ψ, loc_ccv, g, i, j + 0, k + 0))), + ((lerp(ψ, loc_ccv, g, i - 1, j - 1, k + 1), lerp(ψ, loc_ccv, g, i, j - 1, k + 1)), + (lerp(ψ, loc_ccv, g, i - 1, j + 0, k + 1), lerp(ψ, loc_ccv, g, i, j + 0, k + 1)))) + ω.vvc[i, j, k] = ψ2ω(cell_vvc, Δ(g, loc_vvc, i, j, k)) + + # xz edge + cell_vcv = (((lerp(ψ, loc_cvc, g, i - 1, j + 0, k - 1), lerp(ψ, loc_cvc, g, i, j + 0, k - 1)), + (lerp(ψ, loc_cvc, g, i - 1, j + 1, k - 1), lerp(ψ, loc_cvc, g, i, j + 1, k - 1))), + ((lerp(ψ, loc_cvc, g, i - 1, j + 0, k + 0), lerp(ψ, loc_cvc, g, i, j + 0, k + 0)), + (lerp(ψ, loc_cvc, g, i - 1, j + 1, k + 0), lerp(ψ, loc_cvc, g, i, j + 1, k + 0)))) + ω.vcv[i, j, k] = ψ2ω(cell_vcv, Δ(g, loc_vcv, i, j, k)) + + # yz edge + cell_cvv = (((lerp(ψ, loc_vcc, g, i, j - 1, k - 1), lerp(ψ, loc_vcc, g, i + 1, j - 1, k - 1)), + (lerp(ψ, loc_vcc, g, i, j + 0, k - 1), lerp(ψ, loc_vcc, g, i + 1, j + 0, k - 1))), + ((lerp(ψ, loc_vcc, g, i, j - 1, k + 0), lerp(ψ, loc_vcc, g, i + 1, j - 1, k + 0)), + (lerp(ψ, loc_vcc, g, i, j + 0, k + 0), lerp(ψ, loc_vcc, g, i + 1, j + 0, k + 0)))) + ω.cvv[i, j, k] = ψ2ω(cell_cvv, Δ(g, loc_cvv, i, j, k)) +end + +const AnyFieldMask = Union{FieldMask1D, FieldMask2D, FieldMask3D} + +function ω_from_ψ!(arch::Architecture, launch::Launcher, ω::AnyFieldMask, ψ::AbstractField, grid::StructuredGrid) + launch(arch, grid, _ω_from_ψ! => (ω, ψ, grid)) + return +end + +function ω_from_ψ!(arch::DistributedArchitecture, launch::Launcher, ω::AnyFieldMask, ψ::AbstractField, grid::StructuredGrid) + launch(arch, grid, _ω_from_ψ! => (ω, ψ, grid); bc=batch(grid; exchange=Tuple(ω))) + # exchange_halo!(arch, grid, Tuple(ω)...) + return +end diff --git a/src/Models/Diffusion/Heat/Heat.jl b/src/Models/Diffusion/Heat/Heat.jl index 967bd79c..b5a0f0dd 100644 --- a/src/Models/Diffusion/Heat/Heat.jl +++ b/src/Models/Diffusion/Heat/Heat.jl @@ -3,75 +3,57 @@ module Heat export BoundaryCondition, Flux, Value export HeatDiffusionModel, advance_iteration!, advance_timestep! -using FastIce.Grids -using FastIce.Fields -using FastIce.BoundaryConditions -using FastIce.Utils +using Chmy.Architectures +using Chmy.Grids +using Chmy.Fields +using Chmy.BoundaryConditions +using Chmy.KernelLaunch +using Chmy.GridOperators include("kernels.jl") include("boundary_conditions.jl") -struct HeatDiffusionModel{Backend,Grid,BC,Physics,IterParams,Fields} - backend::Backend +struct HeatDiffusionModel{Arch,Grid,BC,Physics,IterParams,Temperature,HeatFlux,KL} + arch::Arch grid::Grid boundary_conditions::BC physics::Physics iter_params::IterParams - fields::Fields + temperature::Temperature + heat_flux::HeatFlux + launcher::KL end -function make_fields_diffusion(backend, grid) - return ( - q = ( - x = Field(backend, grid, (Vertex(), Center(), Center()); halo=1), - y = Field(backend, grid, (Center(), Vertex(), Center()); halo=1), - z = Field(backend, grid, (Center(), Center(), Vertex()); halo=1), - ), - ) -end +function HeatDiffusionModel(; arch, grid, boundary_conditions, physics, iter_params, outer_width=nothing) + temperature = Field(arch, grid, Center()) + heat_flux = VectorField(arch, grid) -function HeatDiffusionModel(; backend, grid, boundary_conditions, physics=nothing, iter_params, fields=nothing, init_fields=nothing) - if isnothing(fields) - diffusion_fields = make_fields_diffusion(backend, grid) - fields = diffusion_fields + if isnothing(outer_width) + outer_width = ntuple(_ -> 2, Val(ndims(grid))) end - if !isnothing(init_fields) - fields = merge(fields, init_fields) - end + launcher = KernelLaunch(arch, grid; outer_width) boundary_conditions = make_field_boundary_conditions(boundary_conditions) - return HeatDiffusionModel(backend, grid, boundary_conditions, physics, iter_params, fields) + return HeatDiffusionModel(arch, grid, boundary_conditions, physics, iter_params, temperature, heat_flux, launcher) end -fields(model::HeatDiffusionModel) = model.fields -grid(model::HeatDiffusionModel) = model.grid - -function advance_iteration!(model::HeatDiffusionModel, t, Δt; async = true) +function advance_iteration!(model::HeatDiffusionModel, t, Δt; async=false) (; T, T_o, q) = model.fields (; Δτ) = model.iter_params - λ_ρCp = model.physics.properties - Δ = NamedTuple{(:x, :y, :z)}(spacing(model.grid)) - nx, ny, nz = size(model.grid) - backend = model.backend - - set_bcs!(bcs) = _apply_bcs!(model.backend, model.grid, model.fields, bcs) + λ_ρCp = model.physics.λ_ρCp + arch = model.arch + grid = model.grid + launch = model.launcher + boundary_conditions = model.boundary_conditions # flux - update_q!(backend, 256, (nx+1, ny+1, nz+1))(q, T, λ_ρCp, Δτ, Δ) - set_bcs!(model.boundary_conditions.flux) - # mass balance - update_T!(backend, 256, (nx, ny, nz))(T, T_o, q, Δt, Δτ, Δ) - set_bcs!(model.boundary_conditions.value) - - async || synchronize(backend) - return -end -function advance_timestep!(model::HeatDiffusionModel, t, Δt) - # TODO + launch(arch, grid, update_q! => (q, T, λ_ρCp, Δτ, grid); bc=boundary_conditions.flux) + launch(arch, grid, update_T! => (T, T_o, q, Δt, Δτ, Δ); bc=boundary_conditions.temperature) + async || synchronize(Architectures.get_backend(arch)) return end diff --git a/src/Models/Diffusion/Heat/boundary_conditions.jl b/src/Models/Diffusion/Heat/boundary_conditions.jl index 21704c60..73a80e14 100644 --- a/src/Models/Diffusion/Heat/boundary_conditions.jl +++ b/src/Models/Diffusion/Heat/boundary_conditions.jl @@ -67,17 +67,3 @@ function make_field_boundary_conditions(bcs) return NamedTuple{(:value, :flux)}(zip(field_bcs...)) end - -function _apply_bcs!(backend, grid, fields, bcs) - field_map = (T = fields.T, - qx = fields.q.x , qy = fields.q.y , qz = fields.q.z) - - ntuple(Val(length(bcs))) do D - if !isnothing(bcs[D]) - fs = Tuple( field_map[f] for f in eachindex(bcs[D]) ) - apply_bcs!(Val(D), backend, grid, fs, values(bcs[D])) - end - end - - return -end diff --git a/src/Models/Diffusion/Heat/kernels.jl b/src/Models/Diffusion/Heat/kernels.jl index db1fb7ad..4650effe 100644 --- a/src/Models/Diffusion/Heat/kernels.jl +++ b/src/Models/Diffusion/Heat/kernels.jl @@ -1,7 +1,5 @@ using KernelAbstractions -using FastIce.GridOperators - @kernel function update_q!(q, T, λ_ρCp, Δτ, Δ) I = @index(Global, Cartesian) @inbounds if checkbounds(Bool, q.x, I) @@ -28,4 +26,4 @@ end ΔTΔt = (T[I] - T_o[I]) / Δt T[I] -= (ΔTΔt + ∇q) * Δτ.T end -end \ No newline at end of file +end diff --git a/src/Models/FullStokes/Isothermal/Isothermal.jl b/src/Models/FullStokes/Isothermal/Isothermal.jl index 5c951a61..17861cdd 100644 --- a/src/Models/FullStokes/Isothermal/Isothermal.jl +++ b/src/Models/FullStokes/Isothermal/Isothermal.jl @@ -3,25 +3,23 @@ module Isothermal export BoundaryCondition, Traction, Velocity, Slip export IsothermalFullStokesModel, advance_iteration!, advance_timestep!, compute_residuals! -using FastIce.Architectures using FastIce.Physics -using FastIce.Grids -using FastIce.Fields -using FastIce.BoundaryConditions -using FastIce.Utils -using FastIce.KernelLaunch -using FastIce.Distributed - -using FastIce.GridOperators + +using Chmy +using Chmy.Architectures +using Chmy.Grids +using Chmy.BoundaryConditions +using Chmy.Distributed +using Chmy.Fields +using Chmy.KernelLaunch +using Chmy.GridOperators + using KernelAbstractions -include("kernels_common.jl") include("kernels_2d.jl") include("kernels_3d.jl") -include("boundary_conditions.jl") - -mutable struct IsothermalFullStokesModel{Arch,Grid,Stress,Velocity,Viscosity,Rheology,Residual,BC,HB,Gravity,SolverParams} +mutable struct IsothermalFullStokesModel{Arch,Grid,Stress,Velocity,Viscosity,Rheology,Residual,BC,Gravity,SolverParams,KL} arch::Arch grid::Grid stress::Stress @@ -31,57 +29,33 @@ mutable struct IsothermalFullStokesModel{Arch,Grid,Stress,Velocity,Viscosity,Rhe gravity::Gravity residual::Residual boundary_conditions::BC - hide_boundaries::HB solver_params::SolverParams + launcher::KL end -function StressFields(arch, grid::CartesianGrid{2}) - (Pr=Field(arch, grid, Center(); halo=1), - # deviatoric stress - τ=(xx=Field(arch, grid, Center(); halo=1), - yy=Field(arch, grid, Center(); halo=1), - xy=Field(arch, grid, (Vertex(), Vertex())))) -end - -function VelocityFields(arch, grid::CartesianGrid{2}) - (x=Field(arch, grid, (Vertex(), Center()); halo=1), - y=Field(arch, grid, (Center(), Vertex()); halo=1)) -end - -function ResidualFields(arch, grid::CartesianGrid{2}) - (r_Pr=Field(arch, grid, Center()), - r_V=(x=Field(arch, grid, (Vertex(), Center())), - y=Field(arch, grid, (Center(), Vertex())))) +struct StressField{Pressure,DeviatoricStress} + P::Pressure + τ::DeviatoricStress end +StressField(arch, grid::StructuredGrid) = StressField(Field(arch, grid, Center()), TensorField(arch, grid)) -function StressFields(arch, grid::CartesianGrid{3}) - (Pr=Field(arch, grid, Center(); halo=1), - # deviatoric stress - τ=(xx=Field(arch, grid, Center(); halo=1), - yy=Field(arch, grid, Center(); halo=1), - zz=Field(arch, grid, Center(); halo=1), - xy=Field(arch, grid, (Vertex(), Vertex(), Center())), - xz=Field(arch, grid, (Vertex(), Center(), Vertex())), - yz=Field(arch, grid, (Center(), Vertex(), Vertex())))) +struct VelocityField{Velocity} + V::Velocity end +VelocityField(arch, grid::StructuredGrid) = VelocityField(VectorField(arch, grid)) -function VelocityFields(arch, grid::CartesianGrid{3}) - (x=Field(arch, grid, (Vertex(), Center(), Center()); halo=1), - y=Field(arch, grid, (Center(), Vertex(), Center()); halo=1), - z=Field(arch, grid, (Center(), Center(), Vertex()); halo=1)) +struct ResidualField{VelocityResidual,PressureResidual} + r_V::VelocityResidual + r_P::PressureResidual end +ResidualField(arch, grid::StructuredGrid) = ResidualField(VectorField(arch, grid), Field(arch, grid, Center())) -function ResidualFields(arch, grid::CartesianGrid{3}) - (r_Pr=Field(arch, grid, Center()), - r_V=(x=Field(arch, grid, (Vertex(), Center(), Center())), - y=Field(arch, grid, (Center(), Vertex(), Center())), - z=Field(arch, grid, (Center(), Center(), Vertex())))) +function ViscosityField(arch, grid::StructuredGrid) + return (η=Field(arch, grid, Center()), + η_next=Field(arch, grid, Center())) end -function ViscosityFields(arch, grid::CartesianGrid) - (η = Field(arch, grid, Center(); halo=1), - η_next = Field(arch, grid, Center(); halo=1)) -end +include("boundary_conditions.jl") function IsothermalFullStokesModel(; arch, grid, @@ -90,18 +64,16 @@ function IsothermalFullStokesModel(; arch, rheology, solver_params=(), outer_width=nothing) - stress = StressFields(arch, grid) - velocity = VelocityFields(arch, grid) - viscosity = ViscosityFields(arch, grid) - residual = ResidualFields(arch, grid) - - boundary_conditions = IsothermalFullStokesBoundaryConditions(arch, grid, stress, velocity, viscosity, residual, boundary_conditions) + stress = StressField(arch, grid) + velocity = VelocityField(arch, grid) + viscosity = ViscosityField(arch, grid) + residual = ResidualField(arch, grid) if isnothing(outer_width) outer_width = ntuple(_ -> 2, Val(ndims(grid))) end - hide_boundaries = HideBoundaries(arch, outer_width) + launcher = Launcher(arch, grid; outer_width) return IsothermalFullStokesModel(arch, grid, @@ -112,57 +84,43 @@ function IsothermalFullStokesModel(; arch, gravity, residual, boundary_conditions, - hide_boundaries, - solver_params) + solver_params, + launcher) end function advance_iteration!(model::IsothermalFullStokesModel, t, Δt) - (; Pr, τ) = model.stress - V = model.velocity - (; η, η_next) = model.viscosity - (; Δτ) = model.solver_params - rheology = model.rheology - ρg = model.gravity - bc = model.boundary_conditions - hide_boundaries = model.hide_boundaries - grid = model.grid - - Δ = spacing(model.grid) - - launch!(model.arch, grid, update_σ! => (Pr, τ, V, η, Δτ, Δ, grid); - location=Center(), expand=1, boundary_conditions=bc.stress, hide_boundaries) - - # merge boundary conditions because viscosity is double-buffered - velocity_bc = dim_side_ntuple(Val(ndims(grid))) do D, S - merge_boundary_conditions(bc.velocity[D][S], bc.viscosity.η_next[D][S]) - end - - launch!(model.arch, grid, update_V! => (V, Pr, τ, η, η_next, rheology, ρg, Δτ, Δ, grid); - location=Vertex(), boundary_conditions=velocity_bc, hide_boundaries) + (; P, τ) = model.stress + (; V) = model.velocity + (; η, η_next) = model.viscosity + (; Δτ) = model.solver_params + arch = model.arch + rheology = model.rheology + ρg = model.gravity + grid = model.grid + launch = model.launcher + + bc1 = batch(grid, model.stress, model.boundary_conditions) + bc2 = merge(batch(grid, model.velocity, model.boundary_conditions; exchange=Tuple(V)), + batch(grid, η_next => Neumann(); exchange=η_next)) + + launch(arch, grid, update_σ! => (P, τ, V, η, Δτ, grid); bc=bc1) + launch(arch, grid, update_V! => (V, P, τ, η, η_next, rheology, ρg, Δτ, grid); bc=bc2) # swap double buffers for viscosity model.viscosity = NamedTuple{keys(model.viscosity)}(reverse(values(model.viscosity))) - bc.viscosity = NamedTuple{keys(bc.viscosity)}(reverse(values(bc.viscosity))) return end function compute_residuals!(model::IsothermalFullStokesModel) - (; Pr, τ) = model.stress - V = model.velocity - (; r_Pr, r_V) = model.residual - ρg = model.gravity - boundary_conditions = model.boundary_conditions.residual - - Δ = spacing(model.grid) - - launch!(model.arch, model.grid, compute_residuals! => (r_V, r_Pr, Pr, τ, V, ρg, Δ, model.grid); - location=Vertex(), boundary_conditions) - return -end - -function advance_timestep!(model::IsothermalFullStokesModel, t, Δt) - # TODO - + (; P, τ) = model.stress + (; V) = model.velocity + (; r_P, r_V) = model.residual + grid = model.grid + ρg = model.gravity + launch = model.launcher + + bc = batch(grid, model.residual, model.boundary_conditions) + launch(model.arch, grid, compute_residuals! => (r_V, r_P, P, τ, V, ρg, grid); bc) return end diff --git a/src/Models/FullStokes/Isothermal/boundary_conditions.jl b/src/Models/FullStokes/Isothermal/boundary_conditions.jl index 285f86df..7f68b693 100644 --- a/src/Models/FullStokes/Isothermal/boundary_conditions.jl +++ b/src/Models/FullStokes/Isothermal/boundary_conditions.jl @@ -12,127 +12,117 @@ end BoundaryCondition{Kind}(components...) where {Kind} = BoundaryCondition{Kind}(components) -function instantiate_boundary_conditions(coords::NTuple{ND}) where {ND} - for (dim, val) in coords - td = remove_dim(Val(val), coords) - TN = ntuple(Val(length(td))) do I - dim2, val2 = td[I] - val2 < val ? Symbol(:τ, dim2, dim) : Symbol(:τ, dim, dim2) - end - N = Symbol(:τ, dim, dim) - - ex1 = Expr(:(=), :Pr, :(DirichletBC{HalfCell}(bc.components[$val]))) - ex2 = Expr(:(=), N, :(DirichletBC{HalfCell}(convert(eltype(bc.components[$val]), 0)))) - ex3 = ntuple(Val(length(TN))) do I - Expr(:(=), TN[I], :(DirichletBC{FullCell}(bc.components[$(td[I][2])]))) - end - - ex_tr = Expr(:tuple, ex1, ex2, ex3...) - - ex_vel = ntuple(length(coords)) do I - kind = I == val ? :(FullCell) : :(HalfCell) - Expr(:(=), Symbol(:V, coords[I][1]), :(DirichletBC{$kind}(bc.components[$I]))) - end - - ex_vel = Expr(:tuple, ex_vel...) - - ex_slip_tr = Expr(:tuple, ex3...) - - ex_slip_vel = Expr(:(=), Symbol(:V, dim), :(DirichletBC{FullCell}(bc.components[$val]))) - ex_slip_vel = Expr(:tuple, ex_slip_vel) - - @eval begin - extract_stress_bc(::Val{$val}, bc::BoundaryCondition{Traction,$ND}) = $ex_tr - extract_stress_bc(::Val{$val}, bc::BoundaryCondition{Velocity,$ND}) = () - extract_stress_bc(::Val{$val}, bc::BoundaryCondition{Slip,$ND}) = $ex_slip_tr - - extract_velocity_bc(::Val{$val}, bc::BoundaryCondition{Traction,$ND}) = () - extract_velocity_bc(::Val{$val}, bc::BoundaryCondition{Velocity,$ND}) = $ex_vel - extract_velocity_bc(::Val{$val}, bc::BoundaryCondition{Slip,$ND}) = $ex_slip_vel - end - end +# 1D boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{1}, stress::StressField, cfd_bc::BoundaryCondition{Traction,1}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[1]) + fields = (stress.P, stress.τ.xx) + bcs = (bc_n, bc_n) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +# 2D boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{2}, stress::StressField, cfd_bc::BoundaryCondition{Traction,2}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[1]) + bc_t = Dirichlet(cfd_bc.components[2]) + fields = (stress.P, stress.τ.xx, stress.τ.xy) + bcs = (bc_n, bc_n, bc_t) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{2}, grid::StructuredGrid{2}, stress::StressField, cfd_bc::BoundaryCondition{Traction,2}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[2]) + bc_t = Dirichlet(cfd_bc.components[1]) + fields = (stress.P, stress.τ.yy, stress.τ.xy) + bcs = (bc_n, bc_n, bc_t) + return batch(side, dim, grid, fields, bcs; kwargs...) end -instantiate_boundary_conditions(((:x, 1), (:y, 2))) -instantiate_boundary_conditions(((:x, 1), (:y, 2), (:z, 3))) +function BoundaryConditions.batch(side, dim::Dim{D}, grid::StructuredGrid{2}, stress::StressField, cfd_bc::BoundaryCondition{Slip,2}; kwargs...) where {D} + bc_t = Dirichlet(cfd_bc.components[3-D]) + fields = tuple(stress.τ.xy) + bcs = tuple(bc_t) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +# 3D boundary conditions -make_batch(::Tuple{}, fields) = nothing +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Traction,3}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[1]) + bc_t1 = Dirichlet(cfd_bc.components[2]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.P, stress.τ.xx, stress.τ.xy, stress.τ.xz) + bcs = (bc_n, bc_n, bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end -function make_batch(bcs::NamedTuple, fields::NamedTuple) - batch_fields = Tuple(fields[name] for name in eachindex(bcs)) - return BoundaryConditionsBatch(batch_fields, values(bcs)) +function BoundaryConditions.batch(side, dim::Dim{2}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Traction,3}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[2]) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.P, stress.τ.yy, stress.τ.xy, stress.τ.yz) + bcs = (bc_n, bc_n, bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) end -@inline function dim_side_ntuple(f::F, ::Val{N}) where {F,N} - ntuple(D -> ntuple(S -> f(D, S), Val(2)), Val(N)) +function BoundaryConditions.batch(side, dim::Dim{3}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Traction,3}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[3]) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[2]) + fields = (stress.P, stress.τ.zz, stress.τ.xz, stress.τ.yz) + bcs = (bc_n, bc_n, bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) end -function StressBoundaryConditions(arch::Architecture{Kind}, ::CartesianGrid{N}, fields, bc) where {Kind,N} - ordering = (:x, :y, :z) - dim_side_ntuple(Val(N)) do D, S - if (Kind == Distributed.DistributedMPI) && has_neighbor(details(arch), D, S) - nothing - else - new_bc = extract_stress_bc(Val(D), bc[ordering[D]][S]) - make_batch(new_bc, fields) - end - end +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Slip,3}; kwargs...) + bc_t1 = Dirichlet(cfd_bc.components[2]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.τ.xy, stress.τ.xz) + bcs = (bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) end -function VelocityBoundaryConditions(arch::Architecture{Kind}, ::CartesianGrid{N}, fields::NamedTuple{names}, bc) where {Kind,N,names} - ordering = (:x, :y, :z) - dim_side_ntuple(Val(N)) do D, S - if (Kind == Distributed.DistributedMPI) && has_neighbor(details(arch), D, S) - new_bc = NamedTuple{names}(ExchangeInfo(Val(S), Val(D), V) for V in fields) - make_batch(new_bc, fields) - else - new_bc = extract_velocity_bc(Val(D), bc[ordering[D]][S]) - make_batch(new_bc, fields) - end - end +function BoundaryConditions.batch(side, dim::Dim{2}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Slip,3}; kwargs...) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.τ.xy, stress.τ.yz) + bcs = (bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) end -function ResidualBoundaryConditions(arch::Architecture{Kind}, ::CartesianGrid{N}, fields::NamedTuple, bc) where {Kind,N} - ordering = (:x, :y, :z) - dim_side_ntuple(Val(N)) do D, S - if (Kind == Distributed.DistributedMPI) && has_neighbor(details(arch), D, S) - nothing - else - if !(bc[ordering[D]][S] isa BoundaryCondition{Traction}) - Vn = Symbol(:r_V, ordering[D]) - BoundaryConditionsBatch((fields[Vn],), (DirichletBC{FullCell}(0.0),)) - else - nothing - end - end - end +function BoundaryConditions.batch(side, dim::Dim{3}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Slip,3}; kwargs...) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[2]) + fields = (stress.τ.xz, stress.τ.yz) + bcs = (bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) end -function ViscosityBoundaryConditions(arch::Architecture{Kind}, ::CartesianGrid{N}, η) where {Kind,N} - dim_side_ntuple(Val(N)) do D, S - if (Kind == Distributed.DistributedMPI) && has_neighbor(details(arch), D, S) - BoundaryConditionsBatch((η,), (ExchangeInfo(Val(S), Val(D), η),)) - else - BoundaryConditionsBatch((η,), (ExtrapolateBC(),)) - end - end +# nD boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{D}, grid::StructuredGrid{N}, velocity::VelocityField, cfd_bc::BoundaryCondition{Velocity,N}; kwargs...) where {D,N} + fields = Tuple(velocity.V) + bcs = ntuple(dim -> Dirichlet(cfd_bc.components[dim]), Val(N)) + return batch(side, dim, grid, fields, bcs; kwargs...) end -mutable struct IsothermalFullStokesBoundaryConditions{Stress,Velocity,Viscosity,Residual} - stress::Stress - velocity::Velocity - viscosity::Viscosity - residual::Residual +function BoundaryConditions.batch(side, dim::Dim{D}, grid::StructuredGrid{N}, velocity::VelocityField, cfd_bc::BoundaryCondition{Slip,N}; kwargs...) where {D,N} + bc_n = Dirichlet(cfd_bc.components[D]) + fields = tuple(velocity.V[D]) + bcs = tuple(bc_n) + return batch(side, dim, grid, fields, bcs; kwargs...) end -function IsothermalFullStokesBoundaryConditions(arch::Architecture, grid::CartesianGrid, stress, velocity, viscosity, residual, bc) - stress_fields = (; Pr=stress.Pr, NamedTuple{Symbol.(:τ, keys(stress.τ))}(values(stress.τ))...) - velocity_fields = NamedTuple{Symbol.(:V, keys(velocity))}(values(velocity)) - residual_fields = NamedTuple{Symbol.(:r_V, keys(residual.r_V))}(values(residual.r_V)) +BoundaryConditions.batch(side, dim, grid, ::StressField, ::BoundaryCondition{Velocity}; kwargs...) = EmptyBatch() +BoundaryConditions.batch(side, dim, grid, ::VelocityField, ::BoundaryCondition{Traction}; kwargs...) = EmptyBatch() + +function BoundaryConditions.batch(side, dim::Dim, grid, residual::ResidualField, ::BoundaryCondition{Traction}; kwargs...) + return EmptyBatch() +end - return IsothermalFullStokesBoundaryConditions(StressBoundaryConditions(arch, grid, stress_fields, bc), - VelocityBoundaryConditions(arch, grid, velocity_fields, bc), - (η = ViscosityBoundaryConditions(arch, grid, viscosity.η), - η_next = ViscosityBoundaryConditions(arch, grid, viscosity.η_next)), - ResidualBoundaryConditions(arch, grid, residual_fields, bc)) +function BoundaryConditions.batch(side, dim::Dim{D}, grid, residual::ResidualField, ::BoundaryCondition; kwargs...) where {D} + fields = tuple(residual.r_V[D]) + bcs = tuple(Dirichlet()) + return batch(side, dim, grid, fields, bcs; kwargs...) end diff --git a/src/Models/FullStokes/Isothermal/kernels_2d.jl b/src/Models/FullStokes/Isothermal/kernels_2d.jl index f13a8e24..58e40459 100644 --- a/src/Models/FullStokes/Isothermal/kernels_2d.jl +++ b/src/Models/FullStokes/Isothermal/kernels_2d.jl @@ -1,75 +1,54 @@ # pseudo-transient update rules -@kernel inbounds=true function update_σ!(Pr, τ, V, η, Δτ, Δ, grid::CartesianGrid{2}, offset=nothing) +@kernel inbounds = true function update_σ!(Pr, τ, V, η, Δτ, g::StructuredGrid{2}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if checkbounds(Bool, Pr, I) - ε̇xx = ∂ᶜx(V.x, I) / Δ.x - ε̇yy = ∂ᶜy(V.y, I) / Δ.y - ∇V = ε̇xx + ε̇yy - # hydrostatic - Pr[I] -= ∇V * η[I] * Δτ.Pr - # deviatoric diagonal - τ.xx[I] -= (τ.xx[I] - 2.0 * η[I] * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx - τ.yy[I] -= (τ.yy[I] - 2.0 * η[I] * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy - end - if checkbounds(Bool, τ.xy, I) - τ.xy[I] -= (τ.xy[I] - avᵛxy(η, I) * (∂ᵛx(V.y, I) / Δ.x + ∂ᵛy(V.x, I) / Δ.y)) * Δτ.τ.xy - end + I += O + # strain rates + ε̇xx = ∂x(V.x, g, I) + ε̇yy = ∂y(V.y, g, I) + ε̇xy = 0.5 * (∂x(V.y, g, I) + ∂y(V.x, g, I)) + # velocity divergence + ∇V = ε̇xx + ε̇yy + # hydrostatic stress + Pr[I] -= ∇V * lerp(η, location(Pr), g, I) * Δτ.Pr + # deviatoric stress + τ.xx[I] -= (τ.xx[I] - 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx + τ.yy[I] -= (τ.yy[I] - 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy + τ.xy[I] -= (τ.xy[I] - 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy) * Δτ.τ.xy end -@kernel inbounds=true function update_V!(V, Pr, τ, η, η_next, rheology, ρg, Δτ, Δ, grid::CartesianGrid{2}, offset=nothing) +@kernel inbounds = true function update_V!(V, Pr, τ, η, η_next, rheology, ρg, Δτ, g::StructuredGrid{2}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if within(grid, V.x, I) - ∂σxx_∂x = (-∂ᵛx(Pr, I) + ∂ᵛx(τ.xx, I)) / Δ.x - ∂τxy_∂y = ∂ᶜy(τ.xy, I) / Δ.y - V.x[I] += (∂σxx_∂x + ∂τxy_∂y - ρg.x[I]) / maxlᵛx(η, I) * Δτ.V.x - end - if within(grid, V.y, I) - ∂σyy_∂y = (-∂ᵛy(Pr, I) + ∂ᵛy(τ.yy, I)) / Δ.y - ∂τxy_∂x = ∂ᶜx(τ.xy, I) / Δ.x - V.y[I] += (∂σyy_∂y + ∂τxy_∂x - ρg.y[I]) / maxlᵛy(η, I) * Δτ.V.y - end - if within(grid, η_next, I) - τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2) + avᶜxy(τ.xy, I)^2) - η_next[I] = rheology(τII, I) - end + I += O + # velocity + V.x[I] += (-∂x(Pr, g, I) + ∂x(τ.xx, g, I) + ∂y(τ.xy, g, I) - ρg.x[I]) / lerp(η, location(V.x), I) * Δτ.V.x + V.y[I] += (-∂y(Pr, g, I) + ∂y(τ.yy, g, I) + ∂x(τ.xy, g, I) - ρg.y[I]) / lerp(η, location(V.x), I) * Δτ.V.y + # rheology + τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2) + lerp(τ.xy, location(η), g, I)^2) + η_next[I] = rheology(τII, I) end # helper kernels -@kernel inbounds=true function compute_τ!(τ, V, η, Δ, grid::CartesianGrid{2}, offset=nothing) +@kernel inbounds = true function compute_τ!(τ, V, η, g::StructuredGrid{2}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if checkbounds(Bool, Pr, I) - ε̇xx = ∂ᶜx(V.x, I) / Δ.x - ε̇yy = ∂ᶜy(V.y, I) / Δ.y - ∇V = ε̇xx + ε̇yy - # deviatoric diagonal - τ.xx[I] = 2.0 * η[I] * (ε̇xx - ∇V / 3.0) - τ.yy[I] = 2.0 * η[I] * (ε̇yy - ∇V / 3.0) - end - if checkbounds(Bool, τ.xy, I) - ε̇xy = 0.5 * (∂ᵛx(V.y, I) / Δ.x + ∂ᵛy(V.x, I) / Δ.y) - τ.xy[I] = 2.0 * avᵛxy(η, I) * ε̇xy - end + I += O + ε̇xx = ∂x(V.x, g, I) + ε̇yy = ∂y(V.y, g, I) + ε̇xy = 0.5 * (∂x(V.y, g, I) + ∂y(V.x, g, I)) + ∇V = ε̇xx + ε̇yy + # deviatoric diagonal + τ.xx[I] = 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0) + τ.yy[I] = 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0) + τ.xy[I] = 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy end -@kernel inbounds=true function compute_residuals!(r_V, r_Pr, Pr, τ, V, ρg, Δ, grid::CartesianGrid{2}, offset=nothing) +@kernel inbounds = true function compute_residuals!(r_V, r_Pr, Pr, τ, V, ρg, g::StructuredGrid{2}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if within(grid, r_Pr, I) - r_Pr[I] = ∂ᶜx(V.x, I) / Δ.x + ∂ᶜy(V.y, I) / Δ.y - end - if within(grid, r_V.x, I) - ∂σxx_∂x = (-∂ᵛx(Pr, I) + ∂ᵛx(τ.xx, I)) / Δ.x - ∂τxy_∂y = ∂ᶜy(τ.xy, I) / Δ.y - r_V.x[I] = ∂σxx_∂x + ∂τxy_∂y - ρg.x[I] - end - if within(grid, r_V.y, I) - ∂σyy_∂y = (-∂ᵛy(Pr, I) + ∂ᵛy(τ.yy, I)) / Δ.y - ∂τxy_∂x = ∂ᶜx(τ.xy, I) / Δ.x - r_V.y[I] = ∂σyy_∂y + ∂τxy_∂x - ρg.y[I] - end + I += O + # pressure + r_Pr[I] = ∂x(V.x, g, I) + ∂y(V.y, g, I) + # velocity + r_V.x[I] = -∂x(Pr, g, I) + ∂x(τ.xx, g, I) + ∂y(τ.xy, g, I) - ρg.x[I] + r_V.y[I] = -∂y(Pr, g, I) + ∂y(τ.yy, g, I) + ∂x(τ.xy, g, I) - ρg.y[I] end diff --git a/src/Models/FullStokes/Isothermal/kernels_3d.jl b/src/Models/FullStokes/Isothermal/kernels_3d.jl index 323bac19..c55e6a73 100644 --- a/src/Models/FullStokes/Isothermal/kernels_3d.jl +++ b/src/Models/FullStokes/Isothermal/kernels_3d.jl @@ -1,109 +1,75 @@ # pseudo-transient update rules -@kernel inbounds = true function update_σ!(Pr, τ, V, η, Δτ, Δ, grid::CartesianGrid{3}, offset=nothing) +@kernel inbounds = true function update_σ!(Pr, τ, V, η, Δτ, g::StructuredGrid{3}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if checkbounds(Bool, Pr, I) - ε̇xx = ∂ᶜx(V.x, I) / Δ.x - ε̇yy = ∂ᶜy(V.y, I) / Δ.y - ε̇zz = ∂ᶜz(V.z, I) / Δ.z - ∇V = ε̇xx + ε̇yy + ε̇zz - # hydrostatic - Pr[I] -= ∇V * η[I] * Δτ.Pr - # deviatoric diagonal - τ.xx[I] -= (τ.xx[I] - 2.0 * η[I] * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx - τ.yy[I] -= (τ.yy[I] - 2.0 * η[I] * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy - τ.zz[I] -= (τ.zz[I] - 2.0 * η[I] * (ε̇zz - ∇V / 3.0)) * Δτ.τ.zz - end + I += O + # strain rates + ε̇xx = ∂x(V.x, g, I) + ε̇yy = ∂y(V.y, g, I) + ε̇zz = ∂z(V.z, g, I) + ε̇xy = 0.5 * (∂x(V.y, g, I) + ∂y(V.x, g, I)) + ε̇xz = 0.5 * (∂x(V.z, g, I) + ∂z(V.x, g, I)) + ε̇yz = 0.5 * (∂y(V.z, g, I) + ∂z(V.y, g, I)) + # velocity divergence + ∇V = ε̇xx + ε̇yy + ε̇zz + # hydrostatic stress + Pr[I] -= ∇V * lerp(η, location(Pr), g, I) * Δτ.Pr + # deviatoric diagonal + τ.xx[I] -= (τ.xx[I] - 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx + τ.yy[I] -= (τ.yy[I] - 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy + τ.zz[I] -= (τ.zz[I] - 2.0 * lerp(η, location(τ.zz), g, I) * (ε̇zz - ∇V / 3.0)) * Δτ.τ.zz # deviatoric off-diagonal - if checkbounds(Bool, τ.xy, I) - τ.xy[I] -= (τ.xy[I] - avᵛxy(η, I) * (∂ᵛx(V.y, I) / Δ.x + ∂ᵛy(V.x, I) / Δ.y)) * Δτ.τ.xy - end - if checkbounds(Bool, τ.xz, I) - τ.xz[I] -= (τ.xz[I] - avᵛxz(η, I) * (∂ᵛx(V.z, I) / Δ.x + ∂ᵛz(V.x, I) / Δ.z)) * Δτ.τ.xz - end - if checkbounds(Bool, τ.yz, I) - τ.yz[I] -= (τ.yz[I] - avᵛyz(η, I) * (∂ᵛy(V.z, I) / Δ.y + ∂ᵛz(V.y, I) / Δ.z)) * Δτ.τ.yz - end + τ.xy[I] -= (τ.xy[I] - 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy) * Δτ.τ.xy + τ.xz[I] -= (τ.xz[I] - 2.0 * lerp(η, location(τ.xz), g, I) * ε̇xz) * Δτ.τ.xz + τ.yz[I] -= (τ.yz[I] - 2.0 * lerp(η, location(τ.yz), g, I) * ε̇yz) * Δτ.τ.yz end -@kernel inbounds = true function update_V!(V, Pr, τ, η, η_next, rheology, ρg, Δτ, Δ, grid::CartesianGrid{3}, offset=nothing) +@kernel inbounds = true function update_V!(V, Pr, τ, η, η_next, rheology, ρg, Δτ, g::StructuredGrid{3}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if within(grid, V.x, I) - ∂σxx_∂x = (-∂ᵛx(Pr, I) + ∂ᵛx(τ.xx, I)) / Δ.x - ∂τxy_∂y = ∂ᶜy(τ.xy, I) / Δ.y - ∂τxz_∂z = ∂ᶜz(τ.xz, I) / Δ.z - V.x[I] += (∂σxx_∂x + ∂τxy_∂y + ∂τxz_∂z - ρg.x[I]) / maxlᵛx(η, I) * Δτ.V.x - end - if within(grid, V.y, I) - ∂σyy_∂y = (-∂ᵛy(Pr, I) + ∂ᵛy(τ.yy, I)) / Δ.y - ∂τxy_∂x = ∂ᶜx(τ.xy, I) / Δ.x - ∂τyz_∂z = ∂ᶜz(τ.yz, I) / Δ.z - V.y[I] += (∂σyy_∂y + ∂τxy_∂x + ∂τyz_∂z - ρg.y[I]) / maxlᵛy(η, I) * Δτ.V.y - end - if within(grid, V.z, I) - ∂σzz_∂z = (-∂ᵛz(Pr, I) + ∂ᵛz(τ.zz, I)) / Δ.z - ∂τxz_∂x = ∂ᶜx(τ.xz, I) / Δ.x - ∂τyz_∂y = ∂ᶜy(τ.yz, I) / Δ.y - V.z[I] += (∂σzz_∂z + ∂τxz_∂x + ∂τyz_∂y - ρg.z[I]) / maxlᵛz(η, I) * Δτ.V.z - end - # update viscosity - if within(grid, η_next, I) - τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2 + τ.zz[I]^2) + - avᶜxy(τ.xy, I)^2 + avᶜxz(τ.xz, I)^2 + avᶜyz(τ.yz, I)^2) - η_next[I] = rheology(τII, I) - end + I += O + # velocity + V.x[I] += (-∂x(Pr, g, I) + ∂x(τ.xx, g, I) + ∂y(τ.xy, g, I) + ∂z(τ.xz, g, I) - ρg.x[I]) / lerp(η, location(V.x), g, I) * Δτ.V.x + V.y[I] += (-∂y(Pr, g, I) + ∂y(τ.yy, g, I) + ∂x(τ.xy, g, I) + ∂z(τ.yz, g, I) - ρg.y[I]) / lerp(η, location(V.y), g, I) * Δτ.V.y + V.z[I] += (-∂z(Pr, g, I) + ∂z(τ.zz, g, I) + ∂x(τ.xz, g, I) + ∂y(τ.yz, g, I) - ρg.z[I]) / lerp(η, location(V.z), g, I) * Δτ.V.z + # rheology + τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2 + τ.zz[I]^2) + + lerp(τ.xy, location(η), g, I)^2 + + lerp(τ.xz, location(η), g, I)^2 + + lerp(τ.yz, location(η), g, I)^2) + η_next[I] = rheology(τII, I) end # helper kernels -@kernel inbounds = true function compute_τ!(τ, V, η, Δ, grid::CartesianGrid{3}, offset=nothing) +@kernel inbounds = true function compute_τ!(τ, V, η, g::StructuredGrid{3}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if checkbounds(Bool, Pr, I) - ε̇xx = ∂ᶜx(V.x, I) / Δ.x - ε̇yy = ∂ᶜy(V.y, I) / Δ.y - ε̇zz = ∂ᶜz(V.z, I) / Δ.z - ∇V = ε̇xx + ε̇yy + ε̇zz - # deviatoric diagonal - τ.xx[I] = 2.0 * η[I] * (ε̇xx - ∇V / 3.0) - τ.yy[I] = 2.0 * η[I] * (ε̇yy - ∇V / 3.0) - τ.zz[I] = 2.0 * η[I] * (ε̇zz - ∇V / 3.0) - end - if checkbounds(Bool, τ.xy, I) - τ.xy[I] = avᵛxy(η, I) * (∂ᵛx(V.y, I) / Δ.x + ∂ᵛy(V.x, I) / Δ.y) - end - if checkbounds(Bool, τ.xz, I) - τ.xz[I] = avᵛxz(η, I) * (∂ᵛx(V.z, I) / Δ.x + ∂ᵛz(V.x, I) / Δ.z) - end - if checkbounds(Bool, τ.yz, I) - τ.yz[I] = avᵛyz(η, I) * (∂ᵛy(V.z, I) / Δ.y + ∂ᵛz(V.y, I) / Δ.z) - end + I += O + # strain rates + ε̇xx = ∂x(V.x, g, I) + ε̇yy = ∂y(V.y, g, I) + ε̇zz = ∂z(V.z, g, I) + ε̇xy = 0.5 * (∂x(V.y, g, I) + ∂y(V.x, g, I)) + ε̇xz = 0.5 * (∂x(V.z, g, I) + ∂z(V.x, g, I)) + ε̇yz = 0.5 * (∂y(V.z, g, I) + ∂z(V.y, g, I)) + # velocity divergence + ∇V = ε̇xx + ε̇yy + ε̇zz + # deviatoric diagonal + τ.xx[I] = 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0) + τ.yy[I] = 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0) + τ.zz[I] = 2.0 * lerp(η, location(τ.zz), g, I) * (ε̇zz - ∇V / 3.0) + # deviatoric off-diagonal + τ.xy[I] = 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy + τ.xz[I] = 2.0 * lerp(η, location(τ.xz), g, I) * ε̇xz + τ.yz[I] = 2.0 * lerp(η, location(τ.yz), g, I) * ε̇yz end -@kernel inbounds = true function compute_residuals!(r_V, r_Pr, Pr, τ, V, ρg, Δ, grid::CartesianGrid{3}, offset=nothing) +@kernel inbounds = true function compute_residuals!(r_V, r_Pr, Pr, τ, V, ρg, g::StructuredGrid{3}, O=Offset()) I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - if within(grid, r_Pr, I) - r_Pr[I] = ∂ᶜx(V.x, I) / Δ.x + ∂ᶜy(V.y, I) / Δ.y + ∂ᶜz(V.z, I) / Δ.z - end - if within(grid, r_V.x, I) - ∂σxx_∂x = (-∂ᵛx(Pr, I) + ∂ᵛx(τ.xx, I)) / Δ.x - ∂τxy_∂y = ∂ᶜy(τ.xy, I) / Δ.y - ∂τxz_∂z = ∂ᶜz(τ.xz, I) / Δ.z - r_V.x[I] = ∂σxx_∂x + ∂τxy_∂y + ∂τxz_∂z - ρg.x[I] - end - if within(grid, r_V.y, I) - ∂σyy_∂y = (-∂ᵛy(Pr, I) + ∂ᵛy(τ.yy, I)) / Δ.y - ∂τxy_∂x = ∂ᶜx(τ.xy, I) / Δ.x - ∂τyz_∂z = ∂ᶜz(τ.yz, I) / Δ.z - r_V.y[I] = ∂σyy_∂y + ∂τxy_∂x + ∂τyz_∂z - ρg.y[I] - end - if within(grid, r_V.z, I) - ∂σzz_∂z = (-∂ᵛz(Pr, I) + ∂ᵛz(τ.zz, I)) / Δ.z - ∂τxz_∂x = ∂ᶜx(τ.xz, I) / Δ.x - ∂τyz_∂y = ∂ᶜy(τ.yz, I) / Δ.y - r_V.z[I] = ∂σzz_∂z + ∂τxz_∂x + ∂τyz_∂y - ρg.z[I] - end + I += O + # pressure + r_Pr[I] = ∂x(V.x, g, I) + ∂y(V.y, g, I) + ∂z(V.z, g, I) + # velocity + r_V.x[I] = -∂x(Pr, g, I) + ∂x(τ.xx, g, I) + ∂y(τ.xy, g, I) + ∂z(τ.xz, g, I) - ρg.x[I] + r_V.y[I] = -∂y(Pr, g, I) + ∂y(τ.yy, g, I) + ∂x(τ.xy, g, I) + ∂z(τ.yz, g, I) - ρg.y[I] + r_V.z[I] = -∂z(Pr, g, I) + ∂z(τ.zz, g, I) + ∂x(τ.xz, g, I) + ∂y(τ.yz, g, I) - ρg.z[I] end diff --git a/src/Models/FullStokes/Isothermal/kernels_common.jl b/src/Models/FullStokes/Isothermal/kernels_common.jl deleted file mode 100644 index 4c38c29a..00000000 --- a/src/Models/FullStokes/Isothermal/kernels_common.jl +++ /dev/null @@ -1,15 +0,0 @@ -Base.@propagate_inbounds @generated function within(grid::CartesianGrid{N}, f::Field{T,N}, I::CartesianIndex{N}) where {T,N} - quote - Base.Cartesian.@nall $N i->I[i] <= size(grid, location(f, Val(i)), i) - end -end - -"Update viscosity using relaxation in log-space. Improves stability of iterative methods" -@kernel function update_η!(η, η_rh, χ, fields, grid, offset=nothing) - I = @index(Global, Cartesian) - isnothing(offset) || (I += offset) - @inbounds begin - ηt = η_rh(grid, I, fields) - η[I] = exp(log(η[I]) * (1 - χ) + log(ηt) * χ) - end -end diff --git a/src/Models/ImmersedBoundaryFullStokes/ImmersedBoundaryFullStokes.jl b/src/Models/ImmersedBoundaryFullStokes/ImmersedBoundaryFullStokes.jl new file mode 100644 index 00000000..eff4f932 --- /dev/null +++ b/src/Models/ImmersedBoundaryFullStokes/ImmersedBoundaryFullStokes.jl @@ -0,0 +1,5 @@ +module ImmersedBoundaryFullStokes + +include("Isothermal/Isothermal.jl") + +end diff --git a/src/Models/ImmersedBoundaryFullStokes/Isothermal/Isothermal.jl b/src/Models/ImmersedBoundaryFullStokes/Isothermal/Isothermal.jl new file mode 100644 index 00000000..653da82c --- /dev/null +++ b/src/Models/ImmersedBoundaryFullStokes/Isothermal/Isothermal.jl @@ -0,0 +1,138 @@ +module Isothermal + +export BoundaryCondition, Traction, Velocity, Slip +export IsothermalImmersedBoundaryFullStokesModel, advance_iteration!, advance_timestep!, compute_residuals! + +using FastIce.Physics + +using Chmy +using Chmy.Architectures +using Chmy.Grids +using Chmy.BoundaryConditions +using Chmy.Distributed +using Chmy.Fields +using Chmy.KernelLaunch +using Chmy.GridOperators + +using KernelAbstractions + +include("null_spaces.jl") +include("kernels_2d.jl") +include("kernels_3d.jl") + +mutable struct IsothermalImmersedBoundaryFullStokesModel{Arch,Grid,Stress,Velocity,Viscosity,Rheology,Residual,BC,LS,FM,Gravity,SolverParams,KL} + arch::Arch + grid::Grid + stress::Stress + velocity::Velocity + viscosity::Viscosity + rheology::Rheology + gravity::Gravity + residual::Residual + boundary_conditions::BC + level_sets::LS + field_masks::FM + solver_params::SolverParams + launcher::KL +end + +struct StressField{Pressure,DeviatoricStress} + P::Pressure + τ::DeviatoricStress +end +StressField(arch, grid::StructuredGrid) = StressField(Field(arch, grid, Center()), TensorField(arch, grid)) + +struct VelocityField{Velocity} + V::Velocity +end +VelocityField(arch, grid::StructuredGrid) = VelocityField(VectorField(arch, grid)) + +struct ResidualField{VelocityResidual,PressureResidual} + r_V::VelocityResidual + r_P::PressureResidual +end +ResidualField(arch, grid::StructuredGrid) = ResidualField(VectorField(arch, grid), Field(arch, grid, Center())) + +function ViscosityField(arch, grid::StructuredGrid) + return (η=Field(arch, grid, Center()), + η_next=Field(arch, grid, Center())) +end + +include("boundary_conditions.jl") + +function IsothermalImmersedBoundaryFullStokesModel(; arch, + grid, + boundary_conditions, + level_sets, + gravity, + rheology, + solver_params=(), + outer_width=nothing) + stress = StressField(arch, grid) + velocity = VelocityField(arch, grid) + viscosity = ViscosityField(arch, grid) + residual = ResidualField(arch, grid) + + field_masks = (na = FieldMask(arch, grid), + ns = FieldMask(arch, grid)) + + if isnothing(outer_width) + outer_width = ntuple(_ -> 4, Val(ndims(grid))) + end + + launcher = Launcher(arch, grid; outer_width) + + return IsothermalImmersedBoundaryFullStokesModel(arch, + grid, + stress, + velocity, + viscosity, + rheology, + gravity, + residual, + boundary_conditions, + level_sets, + field_masks, + solver_params, + launcher) +end + +function advance_iteration!(model::IsothermalImmersedBoundaryFullStokesModel, t, Δt) + (; P, τ) = model.stress + (; V) = model.velocity + (; η, η_next) = model.viscosity + (; Δτ) = model.solver_params + arch = model.arch + rheology = model.rheology + ρg = model.gravity + grid = model.grid + launch = model.launcher + ω = model.field_masks + + bc1 = batch(grid, model.stress, model.boundary_conditions) + bc2 = merge(batch(grid, model.velocity, model.boundary_conditions; exchange=Tuple(V)), + batch(grid, η_next => Neumann(); exchange=η_next)) + + launch(arch, grid, update_σ! => (P, τ, V, η, ω, Δτ, grid); bc=bc1) + launch(arch, grid, update_V! => (V, P, τ, η, η_next, rheology, ρg, ω, Δτ, grid); bc=bc2) + + # swap double buffers for viscosity + model.viscosity = NamedTuple{keys(model.viscosity)}(reverse(values(model.viscosity))) + return +end + +function compute_residuals!(model::IsothermalImmersedBoundaryFullStokesModel) + (; P, τ) = model.stress + (; V) = model.velocity + (; r_P, r_V) = model.residual + grid = model.grid + ρg = model.gravity + launch = model.launcher + ω = model.field_masks + + bc = batch(grid, model.residual, model.boundary_conditions) + launch(model.arch, grid, compute_residuals! => (r_V, r_P, P, τ, V, ρg, ω, grid); bc) + return +end + +end diff --git a/src/Models/ImmersedBoundaryFullStokes/Isothermal/boundary_conditions.jl b/src/Models/ImmersedBoundaryFullStokes/Isothermal/boundary_conditions.jl new file mode 100644 index 00000000..a53d4439 --- /dev/null +++ b/src/Models/ImmersedBoundaryFullStokes/Isothermal/boundary_conditions.jl @@ -0,0 +1,133 @@ +struct Traction end +struct Velocity end +struct Slip end + +""" +Boundary condition for Stokes problem. +""" +struct BoundaryCondition{Kind,N,C} + components::C + BoundaryCondition{Kind}(components::Tuple) where {Kind} = new{Kind,length(components),typeof(components)}(components) +end + +BoundaryCondition{Kind}(components...) where {Kind} = BoundaryCondition{Kind}(components) + +# 1D boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{1}, stress::StressField, cfd_bc::BoundaryCondition{Traction,1}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[1]) + fields = (stress.P, stress.τ.xx) + bcs = (bc_n, bc_n) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +# 2D boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{2}, stress::StressField, cfd_bc::BoundaryCondition{Traction,2}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[1]) + bc_t = Dirichlet(cfd_bc.components[2]) + fields = (stress.P, stress.τ.xx, stress.τ.xy) + bcs = (bc_n, bc_n, bc_t) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{2}, grid::StructuredGrid{2}, stress::StressField, cfd_bc::BoundaryCondition{Traction,2}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[2]) + bc_t = Dirichlet(cfd_bc.components[1]) + fields = (stress.P, stress.τ.yy, stress.τ.xy) + bcs = (bc_n, bc_n, bc_t) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{D}, grid::StructuredGrid{2}, stress::StressField, cfd_bc::BoundaryCondition{Slip,2}; kwargs...) where {D} + bc_t = Dirichlet(cfd_bc.components[3-D]) + fields = tuple(stress.τ.xy) + bcs = tuple(bc_t) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +# 3D boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Traction,3}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[1]) + bc_t1 = Dirichlet(cfd_bc.components[2]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.P, stress.τ.xx, stress.τ.xy, stress.τ.xz) + bcs = (bc_n, bc_n, bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{2}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Traction,3}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[2]) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.P, stress.τ.yy, stress.τ.xy, stress.τ.yz) + bcs = (bc_n, bc_n, bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{3}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Traction,3}; kwargs...) + bc_n = Dirichlet(cfd_bc.components[3]) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[2]) + fields = (stress.P, stress.τ.zz, stress.τ.xz, stress.τ.yz) + bcs = (bc_n, bc_n, bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{1}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Slip,3}; kwargs...) + bc_t1 = Dirichlet(cfd_bc.components[2]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.τ.xy, stress.τ.xz) + bcs = (bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{2}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Slip,3}; kwargs...) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[3]) + fields = (stress.τ.xy, stress.τ.yz) + bcs = (bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{3}, grid::StructuredGrid{3}, stress::StressField, cfd_bc::BoundaryCondition{Slip,3}; kwargs...) + bc_t1 = Dirichlet(cfd_bc.components[1]) + bc_t2 = Dirichlet(cfd_bc.components[2]) + fields = (stress.τ.xz, stress.τ.yz) + bcs = (bc_t1, bc_t2) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +# nD boundary conditions + +function BoundaryConditions.batch(side, dim::Dim{D}, grid::StructuredGrid{N}, velocity::VelocityField, cfd_bc::BoundaryCondition{Velocity,N}; kwargs...) where {D,N} + fields = Tuple(velocity.V) + bcs = ntuple(dim -> Dirichlet(cfd_bc.components[dim]), Val(N)) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{D}, grid::StructuredGrid{N}, velocity::VelocityField, cfd_bc::BoundaryCondition{Slip,N}; kwargs...) where {D,N} + bc_n = Dirichlet(cfd_bc.components[D]) + fields = tuple(velocity.V[D]) + bcs = tuple(bc_n) + return batch(side, dim, grid, fields, bcs; kwargs...) +end + +function BoundaryConditions.batch(side, dim, grid, ::StressField, ::BoundaryCondition{Velocity}; kwargs...) + return batch(side, dim, grid, (), (); kwargs...) +end + +function BoundaryConditions.batch(side, dim, grid, ::VelocityField, ::BoundaryCondition{Traction}; kwargs...) + return batch(side, dim, grid, (), (); kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim, grid, residual::ResidualField, ::BoundaryCondition{Traction}; kwargs...) + return batch(side, dim, grid, (), (); kwargs...) +end + +function BoundaryConditions.batch(side, dim::Dim{D}, grid, residual::ResidualField, ::BoundaryCondition; kwargs...) where {D} + fields = tuple(residual.r_V[D]) + bcs = tuple(Dirichlet()) + return batch(side, dim, grid, fields, bcs; kwargs...) +end diff --git a/src/Models/ImmersedBoundaryFullStokes/Isothermal/kernels_2d.jl b/src/Models/ImmersedBoundaryFullStokes/Isothermal/kernels_2d.jl new file mode 100644 index 00000000..92e30785 --- /dev/null +++ b/src/Models/ImmersedBoundaryFullStokes/Isothermal/kernels_2d.jl @@ -0,0 +1,54 @@ +# pseudo-transient update rules + +@kernel inbounds = true function update_σ!(Pr, τ, V, η, ω, Δτ, g::StructuredGrid{2}, O=Offset()) + I = @index(Global, Cartesian) + I += O + # strain rates + ε̇xx = ∂x(V.x, ω, g, I) + ε̇yy = ∂y(V.y, ω, g, I) + ε̇xy = 0.5 * (∂x(V.y, ω, g, I) + ∂y(V.x, ω, g, I)) + # velocity divergence + ∇V = ε̇xx + ε̇yy + # hydrostatic stress + Pr[I] -= ∇V * lerp(η, location(Pr), g, I) * Δτ.Pr + # deviatoric stress + τ.xx[I] -= (τ.xx[I] - 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx + τ.yy[I] -= (τ.yy[I] - 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy + τ.xy[I] -= (τ.xy[I] - 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy) * Δτ.τ.xy +end + +@kernel inbounds = true function update_V!(V, Pr, τ, η, η_next, rheology, ρg, ω, Δτ, g::StructuredGrid{2}, O=Offset()) + I = @index(Global, Cartesian) + I += O + # velocity + V.x[I] += (-∂x(Pr, ω, g, I) + ∂x(τ.xx, ω, g, I) + ∂y(τ.xy, ω, g, I) - at(ω, location(V.x), I) * ρg.x[I]) / lerp(η, location(V.x), I) * Δτ.V.x + V.y[I] += (-∂y(Pr, ω, g, I) + ∂y(τ.yy, ω, g, I) + ∂x(τ.xy, ω, g, I) - at(ω, location(V.y), I) * ρg.y[I]) / lerp(η, location(V.x), I) * Δτ.V.y + # rheology + τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2) + lerp(τ.xy, location(η), g, I)^2) + η_next[I] = rheology(τII, I) +end + +# helper kernels + +@kernel inbounds = true function compute_τ!(τ, V, η, ω, g::StructuredGrid{2}, O=Offset()) + I = @index(Global, Cartesian) + I += O + ε̇xx = ∂x(V.x, ω, g, I) + ε̇yy = ∂y(V.y, ω, g, I) + ε̇xy = 0.5 * (∂x(V.y, ω, g, I) + ∂y(V.x, ω, g, I)) + ∇V = ε̇xx + ε̇yy + # deviatoric diagonal + τ.xx[I] = 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0) + τ.yy[I] = 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0) + τ.xy[I] = 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy +end + +@kernel inbounds = true function compute_residuals!(r_V, r_Pr, Pr, τ, V, ρg, ω, g::StructuredGrid{2}, O=Offset()) + I = @index(Global, Cartesian) + I += O + # pressure + r_Pr[I] = ∂x(V.x, ω, g, I) + ∂y(V.y, ω, g, I) + # velocity + r_V.x[I] = -∂x(Pr, ω, g, I) + ∂x(τ.xx, ω, g, I) + ∂y(τ.xy, ω, g, I) - at(ω, location(V.x), I) * ρg.x[I] + r_V.y[I] = -∂y(Pr, ω, g, I) + ∂y(τ.yy, ω, g, I) + ∂x(τ.xy, ω, g, I) - at(ω, location(V.y), I) * ρg.y[I] +end diff --git a/src/Models/ImmersedBoundaryFullStokes/Isothermal/kernels_3d.jl b/src/Models/ImmersedBoundaryFullStokes/Isothermal/kernels_3d.jl new file mode 100644 index 00000000..7825a4dd --- /dev/null +++ b/src/Models/ImmersedBoundaryFullStokes/Isothermal/kernels_3d.jl @@ -0,0 +1,143 @@ +# pseudo-transient update rules + +@kernel inbounds = true function update_σ!(Pr, τ, V, η, ω, Δτ, g::StructuredGrid{3}, O=Offset()) + I = @index(Global, Cartesian) + I += O + if isvalid(ω.na, location(Pr), I) && !isempty(ω.ns, location(Pr), I) + # strain rates + ε̇xx = ∂x(V.x, ω.ns, g, I) + ε̇yy = ∂y(V.y, ω.ns, g, I) + ε̇zz = ∂z(V.z, ω.ns, g, I) + # velocity divergence + ∇V = ε̇xx + ε̇yy + ε̇zz + # hydrostatic stress + Pr[I] -= ∇V * lerp(η, location(Pr), g, I) * Δτ.Pr + # deviatoric diagonal + τ.xx[I] -= (τ.xx[I] - 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx + τ.yy[I] -= (τ.yy[I] - 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy + τ.zz[I] -= (τ.zz[I] - 2.0 * lerp(η, location(τ.zz), g, I) * (ε̇zz - ∇V / 3.0)) * Δτ.τ.zz + else + Pr[I] = 0.0 + τ.xx[I] = 0.0 + τ.yy[I] = 0.0 + τ.zz[I] = 0.0 + end + # deviatoric off-diagonal + if isvalid(ω.na, location(τ.xy), I) && !isempty(ω.ns, location(τ.xy), I) + ε̇xy = 0.5 * (∂x(V.y, ω.ns, g, I) + ∂y(V.x, ω.ns, g, I)) + τ.xy[I] -= (τ.xy[I] - 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy) * Δτ.τ.xy + else + τ.xy[I] = 0.0 + end + if isvalid(ω.na, location(τ.xz), I) && !isempty(ω.ns, location(τ.xz), I) + ε̇xz = 0.5 * (∂x(V.z, ω.ns, g, I) + ∂z(V.x, ω.ns, g, I)) + τ.xz[I] -= (τ.xz[I] - 2.0 * lerp(η, location(τ.xz), g, I) * ε̇xz) * Δτ.τ.xz + else + τ.xz[I] = 0.0 + end + if isvalid(ω.na, location(τ.yz), I) && !isempty(ω.ns, location(τ.yz), I) + ε̇yz = 0.5 * (∂y(V.z, ω.ns, g, I) + ∂z(V.y, ω.ns, g, I)) + τ.yz[I] -= (τ.yz[I] - 2.0 * lerp(η, location(τ.yz), g, I) * ε̇yz) * Δτ.τ.yz + else + τ.yz[I] = 0.0 + end +end + +@kernel inbounds = true function update_V!(V, Pr, τ, η, η_next, rheology, ρg, ω, Δτ, g::StructuredGrid{3}, O=Offset()) + I = @index(Global, Cartesian) + I += O + # velocity + if isvalid(ω.ns, location(V.x), I) && !isempty(ω.na, location(V.x), I) + V.x[I] += (-∂x(Pr, ω.na, g, I) + ∂x(τ.xx, ω.na, g, I) + ∂y(τ.xy, ω.na, g, I) + ∂z(τ.xz, ω.na, g, I) - ρg.x[I]) / lerp(η, location(V.x), g, I) / + at(ω.ns, location(V.x), I) * Δτ.V.x + else + V.x[I] = 0.0 + end + if isvalid(ω.ns, location(V.y), I) && !isempty(ω.na, location(V.y), I) + V.y[I] += (-∂y(Pr, ω.na, g, I) + ∂y(τ.yy, ω.na, g, I) + ∂x(τ.xy, ω.na, g, I) + ∂z(τ.yz, ω.na, g, I) - ρg.y[I]) / lerp(η, location(V.y), g, I) / + at(ω.ns, location(V.y), I) * Δτ.V.y + else + V.y[I] = 0.0 + end + if isvalid(ω.ns, location(V.z), I) && !isempty(ω.na, location(V.z), I) + V.z[I] += (-∂z(Pr, ω.na, g, I) + ∂z(τ.zz, ω.na, g, I) + ∂x(τ.xz, ω.na, g, I) + ∂y(τ.yz, ω.na, g, I) - ρg.z[I]) / lerp(η, location(V.z), g, I) / + at(ω.ns, location(V.z), I) * Δτ.V.z + else + V.z[I] = 0.0 + end + # rheology + τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2 + τ.zz[I]^2) + + lerp(τ.xy, location(η), g, I)^2 + + lerp(τ.xz, location(η), g, I)^2 + + lerp(τ.yz, location(η), g, I)^2) + η_next[I] = rheology(τII, I) +end + +# helper kernels + +@kernel inbounds = true function compute_τ!(τ, V, η, ω, g::StructuredGrid{3}, O=Offset()) + I = @index(Global, Cartesian) + I += O + if isvalid(ω.na, location(Pr), I) && !isempty(ω.ns, location(Pr), I) + # strain rates + ε̇xx = ∂x(V.x, ω.ns, g, I) + ε̇yy = ∂y(V.y, ω.ns, g, I) + ε̇zz = ∂z(V.z, ω.ns, g, I) + # velocity divergence + ∇V = ε̇xx + ε̇yy + ε̇zz + # deviatoric diagonal + τ.xx[I] = 2.0 * lerp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0) + τ.yy[I] = 2.0 * lerp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0) + τ.zz[I] = 2.0 * lerp(η, location(τ.zz), g, I) * (ε̇zz - ∇V / 3.0) + else + τ.xx[I] = 0.0 + τ.yy[I] = 0.0 + τ.zz[I] = 0.0 + end + # deviatoric off-diagonal + if isvalid(ω.na, location(τ.xy), I) && !isempty(ω.ns, location(τ.xy), I) + ε̇xy = 0.5 * (∂x(V.y, ω.ns, g, I) + ∂y(V.x, ω.ns, g, I)) + τ.xy[I] = 2.0 * lerp(η, location(τ.xy), g, I) * ε̇xy + else + τ.xy[I] = 0.0 + end + if isvalid(ω.na, location(τ.xz), I) && !isempty(ω.ns, location(τ.xz), I) + ε̇xz = 0.5 * (∂x(V.z, ω.ns, g, I) + ∂z(V.x, ω.ns, g, I)) + τ.xz[I] = 2.0 * lerp(η, location(τ.xz), g, I) * ε̇xz + else + τ.xz[I] = 0.0 + end + if isvalid(ω.na, location(τ.yz), I) && !isempty(ω.ns, location(τ.yz), I) + ε̇yz = 0.5 * (∂y(V.z, ω.ns, g, I) + ∂z(V.y, ω.ns, g, I)) + τ.yz[I] = 2.0 * lerp(η, location(τ.yz), g, I) * ε̇yz + else + τ.yz[I] = 0.0 + end +end + +@kernel inbounds = true function compute_residuals!(r_V, r_Pr, Pr, τ, V, ρg, ω, g::StructuredGrid{3}, O=Offset()) + I = @index(Global, Cartesian) + I += O + # pressure + if isvalid(ω.na, location(Pr), I) && !isempty(ω.ns, location(Pr), I) + r_Pr[I] = ∂x(V.x, ω.ns, g, I) + ∂y(V.y, ω.ns, g, I) + ∂z(V.z, ω.ns, g, I) + else + r_Pr[I] = 0.0 + end + # velocity + if isvalid(ω.ns, location(V.x), I) && !isempty(ω.na, location(V.x), I) + r_V.x[I] = -∂x(Pr, ω.na, g, I) + ∂x(τ.xx, ω.na, g, I) + ∂y(τ.xy, ω.na, g, I) + ∂z(τ.xz, ω.na, g, I) - ρg.x[I] + else + r_V.x[I] = 0.0 + end + if isvalid(ω.ns, location(V.y), I) && !isempty(ω.na, location(V.y), I) + r_V.y[I] = -∂y(Pr, ω.na, g, I) + ∂y(τ.yy, ω.na, g, I) + ∂x(τ.xy, ω.na, g, I) + ∂z(τ.yz, ω.na, g, I) - ρg.y[I] + else + r_V.y[I] = 0.0 + end + if isvalid(ω.ns, location(V.z), I) && !isempty(ω.na, location(V.z), I) + r_V.z[I] = -∂z(Pr, ω.na, g, I) + ∂z(τ.zz, ω.na, g, I) + ∂x(τ.xz, ω.na, g, I) + ∂y(τ.yz, ω.na, g, I) - ρg.z[I] + else + r_V.z[I] = 0.0 + end +end diff --git a/src/Models/ImmersedBoundaryFullStokes/Isothermal/null_spaces.jl b/src/Models/ImmersedBoundaryFullStokes/Isothermal/null_spaces.jl new file mode 100644 index 00000000..3f23b713 --- /dev/null +++ b/src/Models/ImmersedBoundaryFullStokes/Isothermal/null_spaces.jl @@ -0,0 +1,43 @@ +flipped(loc::NTuple{N,Location}, ::Dim{D}) where {N,D} = ntuple(dim -> dim == D ? flip(loc[dim]) : loc[dim], Val(N)) + +Chmy.@add_cartesian function left(ω::AbstractMask{T,N}, loc::NTuple{N,Location}, dim::Dim{D}, I::Vararg{Integer,N}) where {T,N,D} + Il = GridOperators.il(flip(loc[D]), loc[D], dim, I...) + return at(ω, flipped(loc, dim), Il...) +end + +Chmy.@add_cartesian function right(ω::AbstractMask{T,N}, loc::NTuple{N,Location}, dim::Dim{D}, I::Vararg{Integer,N}) where {T,N,D} + Ir = GridOperators.ir(flip(loc[D]), loc[D], dim, I...) + return at(ω, flipped(loc, dim), Ir...) +end + +Base.@propagate_inbounds @generated function generic_isnullspace(ω::AbstractMask{T,N}, loc::NTuple{N,Location}, I::Vararg{Integer,N}) where {T,N} + quote + @inline + Base.Cartesian.@nany $N D -> (left(ω, loc, Dim(D), I...) < 1e-6) || (right(ω, loc, Dim(D), I...) < 1e-6) + end +end + +Base.@propagate_inbounds function isnullspace(ω::AbstractMask{T,3}, loc::NTuple{3,Location}, I::Vararg{Integer,3}) where {T} + return generic_isnullspace(ω, loc, I...) +end + +Base.@propagate_inbounds function isnullspace(ω::AbstractMask{T,3}, loc::Tuple{Vertex,Vertex,Center}, I::Vararg{Integer,3}) where {T} + return (left(ω, loc, Dim(1), I...) < 1e-6) || (right(ω, loc, Dim(1), I...) < 1e-6) || + (left(ω, loc, Dim(2), I...) < 1e-6) || (right(ω, loc, Dim(2), I...) < 1e-6) +end + +Base.@propagate_inbounds function isnullspace(ω::AbstractMask{T,3}, loc::Tuple{Vertex,Center,Vertex}, I::Vararg{Integer,3}) where {T} + return (left(ω, loc, Dim(1), I...) < 1e-6) || (right(ω, loc, Dim(1), I...) < 1e-6) || + (left(ω, loc, Dim(3), I...) < 1e-6) || (right(ω, loc, Dim(3), I...) < 1e-6) +end + +Base.@propagate_inbounds function isnullspace(ω::AbstractMask{T,3}, loc::Tuple{Center,Vertex,Vertex}, I::Vararg{Integer,3}) where {T} + return (left(ω, loc, Dim(2), I...) < 1e-6) || (right(ω, loc, Dim(2), I...) < 1e-6) || + (left(ω, loc, Dim(3), I...) < 1e-6) || (right(ω, loc, Dim(3), I...) < 1e-6) +end + +Base.@propagate_inbounds isnullspace(ω::AbstractMask, loc, I::CartesianIndex) = isnullspace(ω, loc, Tuple(I)...) + +Chmy.@add_cartesian isempty(ω::AbstractMask{T,N}, loc, I::Vararg{Integer,N}) where {T,N} = at(ω, loc, I...) < 1e-6 + +Chmy.@add_cartesian isvalid(ω::AbstractMask{T,N}, loc, I::Vararg{Integer,N}) where {T,N} = !(isnullspace(ω, loc, I...) || isempty(ω, loc, I...)) diff --git a/src/Models/Models.jl b/src/Models/Models.jl index 7668b4b8..3e821ac7 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -2,5 +2,6 @@ module Models include("Diffusion/Diffusion.jl") include("FullStokes/FullStokes.jl") +include("ImmersedBoundaryFullStokes/ImmersedBoundaryFullStokes.jl") end diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl deleted file mode 100644 index e159f4c1..00000000 --- a/src/Utils/Utils.jl +++ /dev/null @@ -1,36 +0,0 @@ -module Utils - -export remove_dim, insert_dim -export split_ndrange -export Worker - -using FastIce.Fields -using FastIce.Architectures -using KernelAbstractions - -include("workers.jl") -include("split_ndrange.jl") - -# Returns a copy of the tuple `A` with element in position `D` removed -@inline remove_dim(::Val{D}, A::NTuple{N}) where {D,N} = ntuple(I -> I < D ? A[I] : A[I+1], Val(N - 1)) -@inline remove_dim(::Val{1}, I::NTuple{1}) = 1 - -# Same as `remove_dim`, but for `CartesianIndex` -@inline remove_dim(dim, I::CartesianIndex) = remove_dim(dim, Tuple(I)) |> CartesianIndex - -# Returns a copy of tuple `A`, but inserts `i` into position `D` -@inline insert_dim(::Val{D}, A::NTuple{N}, i) where {D,N} = - ntuple(Val(N + 1)) do I - if I < D - A[I] - elseif I == D - i - else - A[I-1] - end - end - -# Same as `insert_dim`, but for `CartesianIndex` -@inline insert_dim(dim, A::CartesianIndex, i) = insert_dim(dim, Tuple(A), i) |> CartesianIndex - -end diff --git a/src/Utils/split_ndrange.jl b/src/Utils/split_ndrange.jl deleted file mode 100644 index a50d017b..00000000 --- a/src/Utils/split_ndrange.jl +++ /dev/null @@ -1,26 +0,0 @@ -@inline outer_subrange(nr, bw, D, ::Val{1}) = 1:bw[D] -@inline outer_subrange(nr, bw, D, ::Val{2}) = (size(nr, D)-bw[D]+1):size(nr, D) -@inline inner_subrange(nr, bw, D) = (bw[D]+1):(size(nr, D)-bw[D]) - -function ndsubrange(ndrange::CartesianIndices{N}, ndwidth, I, ::Val{S}) where {N,S} - return ntuple(Val(N)) do D - if D < I - 1:size(ndrange, D) - elseif D == I - outer_subrange(ndrange, ndwidth, D, Val(S)) - else - inner_subrange(ndrange, ndwidth, D) - end - end -end - -@inline split_ndrange(ndrange, ndwidth) = split_ndrange(CartesianIndices(ndrange), ndwidth) - -function split_ndrange(ndrange::CartesianIndices{N}, ndwidth::NTuple{N,<:Integer}) where {N} - @assert all(size(ndrange) .> ndwidth .* 2) - ndinner = ntuple(D -> inner_subrange(ndrange, ndwidth, D), Val(N)) - inner_range = ndrange[ndinner...] - outer_ranges = ntuple(D -> (ndrange[ndsubrange(ndrange, ndwidth, D, Val(1))...], - ndrange[ndsubrange(ndrange, ndwidth, D, Val(2))...]), Val(N)) - return inner_range, outer_ranges -end diff --git a/src/Utils/workers.jl b/src/Utils/workers.jl deleted file mode 100644 index d30e6a35..00000000 --- a/src/Utils/workers.jl +++ /dev/null @@ -1,45 +0,0 @@ -mutable struct Worker - src::Channel - out::Base.Event - task::Task - - function Worker(; setup=nothing, teardown=nothing) - src = Channel() - out = Base.Event(true) - task = Threads.@spawn begin - isnothing(setup) || Base.invokelatest(setup) - try - for work in src - Base.invokelatest(work) - notify(out) - end - finally - isnothing(teardown) || Base.invokelatest(teardown) - end - end - errormonitor(task) - return new(src, out, task) - end -end - -function Base.close(p::Worker) - close(p.src) - wait(p.task) - return -end - -Base.isopen(p::Worker) = isopen(p.src) - -function Base.put!(work::F, p::Worker) where {F} - put!(p.src, work) - return -end - -function Base.wait(p::Worker) - if isopen(p.src) - wait(p.out) - else - error("Worker is not running") - end - return -end diff --git a/src/Writers.jl b/src/Writers.jl index 0a58c7fc..d84f432b 100644 --- a/src/Writers.jl +++ b/src/Writers.jl @@ -2,44 +2,40 @@ module Writers export write_h5, write_xdmf -using FastIce.Architectures -using FastIce.Distributed -using FastIce.Grids -using FastIce.Fields +using Chmy.Architectures +using Chmy.Distributed +using Chmy.Grids +using Chmy.Fields using HDF5 using LightXML using MPI """ - write_h5(arch::Architecture{DistributedMPI}, grid::CartesianGrid, path, fields) + write_h5(arch, grid, path, fields) -Write output `fields` in HDF5 format to a file on `path` for global `grid` on distributed `arch`. +Write output `fields` in HDF5 format to a file on `path` for global `grid`. """ -function write_h5(arch::Architecture{DistributedMPI}, grid::CartesianGrid, path, fields) +function write_h5(::Architecture, grid::UniformGrid, path, fields) + I = CartesianIndices(size(grid, Center())) + h5open(path, "w") do io + write_dset(io, fields, size(grid, Center()), I.indices) + end + return +end + +function write_h5(arch::DistributedArchitecture, grid::UniformGrid, path, fields) HDF5.has_parallel() || @warn("HDF5 has no parallel support.") - topo = details(arch) - comm = cartesian_communicator(topo) - coords = coordinates(topo) - sz = size(local_grid(grid, topo)) + topo = topology(arch) + comm = cart_comm(topo) + coords = cart_coords(topo) + sz = size(grid, Center()) + global_grid_size = dims(topo) .* sz c1 = coords .* sz .+ 1 |> CartesianIndex c2 = (coords .+ 1) .* sz |> CartesianIndex I = c1:c2 h5open(path, "w", comm, MPI.Info()) do io - write_dset(io, fields, size(grid), I.indices) - end - return -end - -""" - write_h5(arch::Architecture, grid::CartesianGrid, path, fields) - -Write output `fields` in HDF5 format to a file on `path`. -""" -function write_h5(arch::Architecture, grid::CartesianGrid, path, fields) - I = CartesianIndices(size(grid)) - h5open(path, "w") do io - write_dset(io, fields, size(grid), I.indices) + write_dset(io, fields, global_grid_size, I.indices) end return end @@ -53,16 +49,15 @@ function write_dset(io, fields, grid_size, inds) end """ - write_xdmf(arch::Architecture{DistributedMPI}, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) + write_xdmf(arch, grid, path, fields, h5_names, timesteps=Float64(0.0)) -Write Xdmf metadata to `path` for corresponding `h5_names` and `fields` for global `grid` on distributed `arch`. +Write Xdmf metadata to `path` for corresponding `h5_names` and `fields` for global `grid`. Saving time-dependant data can be achieved upon passing a vector to `h5_names` and `timesteps`. """ -function write_xdmf(arch::Architecture{DistributedMPI}, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) - topo = details(arch) - grid_size = size(grid) +function write_xdmf(::Architecture, grid::UniformGrid, path, fields, h5_names, timesteps=Float64(0.0)) + grid_size = size(grid, Center()) grid_spacing = spacing(grid) - grid_origin = origin(local_grid(grid, topo)) + grid_origin = origin(grid, Center()) xdoc = generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) @@ -70,18 +65,13 @@ function write_xdmf(arch::Architecture{DistributedMPI}, grid::CartesianGrid, pat return end -""" - write_xdmf(arch::Architecture, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) - -Write Xdmf metadata to `path` for corresponding `h5_names` and `fields`. -Saving time-dependant data can be achieved upon passing a vector to `h5_names` and `timesteps`. -""" -function write_xdmf(arch::Architecture, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) - grid_size = size(grid) +function write_xdmf(arch::DistributedArchitecture, grid::UniformGrid, path, fields, h5_names, timesteps=Float64(0.0)) + topo = topology(arch) + global_grid_size = dims(topo) .* size(grid, Center()) grid_spacing = spacing(grid) - grid_origin = origin(grid) + grid_origin = origin(grid, Center()) - xdoc = generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) + xdoc = generate_xdmf(global_grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) save_file(xdoc, path) return @@ -113,7 +103,7 @@ function generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, t xorig = new_child(xgeom, "DataItem") set_attribute(xorig, "Format", "XML") set_attribute(xorig, "NumberType", "Float") - set_attribute(xorig, "Dimensions", "$(length(grid_size)) ") + set_attribute(xorig, "Dimensions", "$(length(grid_size))") add_text(xorig, join(reverse(grid_origin), ' ')) xdr = new_child(xgeom, "DataItem") @@ -123,6 +113,8 @@ function generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, t add_text(xdr, join(reverse(grid_spacing), ' ')) h5_path = h5_names[it] + @assert endswith(h5_path, ".h5") + for (name, _) in fields create_xdmf_attribute(xgrid, h5_path, name, grid_size) end diff --git a/test/Project.toml b/test/Project.toml index a2881bc9..4b8fdc11 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" @@ -11,5 +12,5 @@ AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" [compat] -AMDGPU = "0.8" +AMDGPU = "0.8, 0.9" CUDA = "5" diff --git a/test/runtests.jl b/test/runtests.jl index 32dd4fe2..11e776c5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,10 +3,10 @@ using FastIce using Pkg -excludedfiles = ["test_excluded.jl"] +excludedfiles = ["test_excluded.jl"] # tmp as WIP # distributed -test_distributed = ["test_distributed_2D.jl", "test_distributed_3D.jl", "test_distributed_writers_3D.jl"] +test_distributed = ["test_distributed_writers_3D.jl"] using MPI nprocs_2D = 4 nprocs_3D = 8 @@ -35,7 +35,6 @@ function parse_flags!(args, flag; default=nothing, typ=typeof(default)) end function runtests() - exename = joinpath(Sys.BINDIR, Base.julia_exename()) testdir = pwd() istest(f) = endswith(f, ".jl") && startswith(basename(f), "test_") testfiles = sort(filter(istest, vcat([joinpath.(root, files) for (root, dirs, files) in walkdir(testdir)]...))) diff --git a/test/test_bc.jl b/test/test_bc.jl deleted file mode 100644 index 8e4bb742..00000000 --- a/test/test_bc.jl +++ /dev/null @@ -1,93 +0,0 @@ -include("common.jl") - -using FastIce.Architectures -using FastIce.BoundaryConditions -using FastIce.Fields -using FastIce.Grids - -for backend in backends - @testset "$(basename(@__FILE__)) (backend: $backend)" begin - @testset "boundary conditions" begin - nx, ny, nz = 2, 2, 2 - grid = CartesianGrid(; origin=(0.0, 0.0, 0.0), extent=(1.0, 1.0, 1.0), size=(nx, ny, nz)) - field = Field(backend, grid, Center(); halo=1) - arch = Architecture(backend) - @testset "value" begin - @testset "x-dim" begin - data(field) .= 0.0 - west_bc = BoundaryConditionsBatch((field,), (DirichletBC{HalfCell}(1.0),)) - east_bc = BoundaryConditionsBatch((field,), (DirichletBC{FullCell}(0.5),)) - apply_boundary_conditions!(Val(1), Val(1), arch, grid, west_bc; async=false) - apply_boundary_conditions!(Val(2), Val(1), arch, grid, east_bc; async=false) - @test all((parent(field)[1, 2:ny+1, 2:nz+1] .+ parent(field)[2, 2:ny+1, 2:nz+1]) ./ 2 .≈ - west_bc.conditions[1].condition) - @test all(parent(field)[nx+2, 2:ny+1, 2:nz+1] .≈ east_bc.conditions[1].condition) - end - @testset "y-dim" begin - data(field) .= 0.0 - south_bc = BoundaryConditionsBatch((field,), (DirichletBC{HalfCell}(1.0),)) - north_bc = BoundaryConditionsBatch((field,), (DirichletBC{FullCell}(0.5),)) - apply_boundary_conditions!(Val(1), Val(2), arch, grid, south_bc; async=false) - apply_boundary_conditions!(Val(2), Val(2), arch, grid, north_bc; async=false) - @test all((parent(field)[2:nx+1, 1, 2:nz+1] .+ parent(field)[2:nx+1, 2, 2:nz+1]) ./ 2 .≈ - south_bc.conditions[1].condition) - @test all(parent(field)[2:nx+1, ny+2, 2:nz+1] .≈ north_bc.conditions[1].condition) - end - @testset "z-dim" begin - data(field) .= 0.0 - bottom_bc = BoundaryConditionsBatch((field,), (DirichletBC{HalfCell}(1.0),)) - top_bc = BoundaryConditionsBatch((field,), (DirichletBC{FullCell}(0.5),)) - apply_boundary_conditions!(Val(1), Val(3), arch, grid, bottom_bc; async=false) - apply_boundary_conditions!(Val(2), Val(3), arch, grid, top_bc; async=false) - @test all((parent(field)[2:nx+1, 2:ny+1, 1] .+ parent(field)[2:nx+1, 2:ny+1, 2]) ./ 2 .≈ - bottom_bc.conditions[1].condition) - @test all(parent(field)[2:nx+1, 2:ny+1, nz+2] .≈ top_bc.conditions[1].condition) - end - end - @testset "array" begin - @testset "x-dim" begin - data(field) .= 0.0 - bc_array_west = KernelAbstractions.allocate(backend, Float64, (size(grid, 2), size(grid, 3))) - bc_array_east = KernelAbstractions.allocate(backend, Float64, (size(grid, 2), size(grid, 3))) - bc_array_west .= 1.0 - bc_array_east .= 0.5 - west_bc = BoundaryConditionsBatch((field,), (DirichletBC{HalfCell}(bc_array_west),)) - east_bc = BoundaryConditionsBatch((field,), (DirichletBC{FullCell}(bc_array_east),)) - apply_boundary_conditions!(Val(1), Val(1), arch, grid, west_bc; async=false) - apply_boundary_conditions!(Val(2), Val(1), arch, grid, east_bc; async=false) - @test all((parent(field)[1, 2:ny+1, 2:nz+1] .+ parent(field)[2, 2:ny+1, 2:nz+1]) ./ 2 .≈ - west_bc.conditions[1].condition) - @test all(parent(field)[nx+2, 2:ny+1, 2:nz+1] .≈ east_bc.conditions[1].condition) - end - @testset "y-dim" begin - data(field) .= 0.0 - bc_array_south = KernelAbstractions.allocate(backend, Float64, (size(grid, 2), size(grid, 3))) - bc_array_north = KernelAbstractions.allocate(backend, Float64, (size(grid, 2), size(grid, 3))) - bc_array_south .= 1.0 - bc_array_north .= 0.5 - south_bc = BoundaryConditionsBatch((field,), (DirichletBC{HalfCell}(bc_array_south),)) - north_bc = BoundaryConditionsBatch((field,), (DirichletBC{FullCell}(bc_array_north),)) - apply_boundary_conditions!(Val(1), Val(2), arch, grid, south_bc; async=false) - apply_boundary_conditions!(Val(2), Val(2), arch, grid, north_bc; async=false) - @test all((parent(field)[2:nx+1, 1, 2:nz+1] .+ parent(field)[2:nx+1, 2, 2:nz+1]) ./ 2 .≈ - south_bc.conditions[1].condition) - @test all(parent(field)[2:nx+1, ny+2, 2:nz+1] .≈ north_bc.conditions[1].condition) - end - @testset "z-dim" begin - data(field) .= 0.0 - bc_array_bot = KernelAbstractions.allocate(backend, Float64, (size(grid, 2), size(grid, 3))) - bc_array_top = KernelAbstractions.allocate(backend, Float64, (size(grid, 2), size(grid, 3))) - bc_array_bot .= 1.0 - bc_array_top .= 0.5 - bottom_bc = BoundaryConditionsBatch((field,), (DirichletBC{HalfCell}(bc_array_bot),)) - top_bc = BoundaryConditionsBatch((field,), (DirichletBC{FullCell}(bc_array_top),)) - apply_boundary_conditions!(Val(1), Val(3), arch, grid, bottom_bc; async=false) - apply_boundary_conditions!(Val(2), Val(3), arch, grid, top_bc; async=false) - @test all((parent(field)[2:nx+1, 2:ny+1, 1] .+ parent(field)[2:nx+1, 2:ny+1, 2]) ./ 2 .≈ - bottom_bc.conditions[1].condition) - @test all(parent(field)[2:nx+1, 2:ny+1, nz+2] .≈ top_bc.conditions[1].condition) - end - end - end - end -end diff --git a/test/test_chmy_grids.jl b/test/test_chmy_grids.jl new file mode 100644 index 00000000..e6c494d9 --- /dev/null +++ b/test/test_chmy_grids.jl @@ -0,0 +1,163 @@ +include("common.jl") + +using Chmy +using Chmy.Architectures +using Chmy.Grids + +@testset "$(basename(@__FILE__)) (backend: CPU)" begin + @testset "common" begin + @test flip(Center()) == Vertex() + @test flip(Vertex()) == Center() + end + + @testset "grids" begin + arch = Arch(CPU()) + nx, ny = 5, 20 + @testset "uniform grids" begin + grid = UniformGrid(arch; origin=(-1, -2), extent=(2, 4), dims=(nx, ny)) + @test grid isa UniformGrid + + # connectivity + @test connectivity(grid, Dim(1), Side(1)) isa Bounded + @test connectivity(grid, Dim(1), Side(2)) isa Bounded + @test connectivity(grid, Dim(2), Side(1)) isa Bounded + @test connectivity(grid, Dim(2), Side(2)) isa Bounded + + # axes + @test axis(grid, Dim(1)) == grid.axes[1] + @test axis(grid, Dim(2)) == grid.axes[2] + + @testset "size" begin + @test size(grid, Center()) == (nx, ny) + @test size(grid, Vertex()) == (nx + 1, ny + 1) + @test size(grid, (Center(), Vertex())) == (nx, ny + 1) + @test size(grid, (Vertex(), Center())) == (nx + 1, ny) + # repeating locations + @test size(grid, Center()) == size(grid, (Center(), Center())) + @test size(grid, Vertex()) == size(grid, (Vertex(), Vertex())) + end + + @testset "bounds" begin + @test all(bounds(grid, Vertex(), Dim(1)) .≈ (-1.0, 1.0)) + @test all(bounds(grid, Vertex(), Dim(2)) .≈ (-2.0, 2.0)) + @test all(bounds(grid, Center(), Dim(1)) .≈ (-0.8, 0.8)) + @test all(bounds(grid, Center(), Dim(2)) .≈ (-1.9, 1.9)) + end + + @testset "extent" begin + # one location + @test extent(grid, Vertex(), Dim(1)) ≈ 2.0 + @test extent(grid, Vertex(), Dim(2)) ≈ 4.0 + @test extent(grid, Center(), Dim(1)) ≈ 1.6 + @test extent(grid, Center(), Dim(2)) ≈ 3.8 + # many locations + @test all(extent(grid, (Vertex(), Vertex())) .≈ (2.0, 4.0)) + @test all(extent(grid, (Center(), Center())) .≈ (1.6, 3.8)) + @test all(extent(grid, (Center(), Vertex())) .≈ (1.6, 4.0)) + @test all(extent(grid, (Vertex(), Center())) .≈ (2.0, 3.8)) + # repeating locations + @test extent(grid, Vertex()) == extent(grid, (Vertex(), Vertex())) + @test extent(grid, Center()) == extent(grid, (Center(), Center())) + end + + @testset "origin" begin + # one location + @test origin(grid, Vertex(), Dim(1)) ≈ -1 + @test origin(grid, Vertex(), Dim(2)) ≈ -2 + @test origin(grid, Center(), Dim(1)) ≈ -0.8 + @test origin(grid, Center(), Dim(2)) ≈ -1.9 + # many locations + @test all(origin(grid, (Vertex(), Vertex())) .≈ (-1.0, -2.0)) + @test all(origin(grid, (Center(), Center())) .≈ (-0.8, -1.9)) + @test all(origin(grid, (Center(), Vertex())) .≈ (-0.8, -2.0)) + @test all(origin(grid, (Vertex(), Center())) .≈ (-1.0, -1.9)) + # repeating locations + @test origin(grid, (Vertex(), Vertex())) == origin(grid, Vertex()) + @test origin(grid, (Center(), Center())) == origin(grid, Center()) + end + + @testset "spacing" begin + @test Δ == spacing + # one location + @test spacing(grid, Vertex(), Dim(1), 1) ≈ 0.4 + @test spacing(grid, Vertex(), Dim(2), 1) ≈ 0.2 + @test spacing(grid, Center(), Dim(1), 1) ≈ 0.4 + @test spacing(grid, Center(), Dim(2), 1) ≈ 0.2 + # many locations + @test all(spacing(grid, (Center(), Center()), 1, 1) .≈ (0.4, 0.2)) + @test all(spacing(grid, (Vertex(), Vertex()), 1, 1) .≈ (0.4, 0.2)) + @test all(spacing(grid, (Center(), Vertex()), 1, 1) .≈ (0.4, 0.2)) + @test all(spacing(grid, (Vertex(), Center()), 1, 1) .≈ (0.4, 0.2)) + # repeating locations + @test spacing(grid, Vertex(), 1, 1) == spacing(grid, (Vertex(), Vertex()), 1, 1) + @test spacing(grid, Center(), 1, 1) == spacing(grid, (Center(), Center()), 1, 1) + end + + @testset "inverse spacing" begin + @test iΔ == inv_spacing + # one location + @test inv_spacing(grid, Vertex(), Dim(1), 1) ≈ 2.5 + @test inv_spacing(grid, Vertex(), Dim(2), 1) ≈ 5.0 + @test inv_spacing(grid, Center(), Dim(1), 1) ≈ 2.5 + @test inv_spacing(grid, Center(), Dim(2), 1) ≈ 5.0 + # many locations + @test all(inv_spacing(grid, (Center(), Center()), 1, 1) .≈ (2.5, 5.0)) + @test all(inv_spacing(grid, (Vertex(), Vertex()), 1, 1) .≈ (2.5, 5.0)) + @test all(inv_spacing(grid, (Center(), Vertex()), 1, 1) .≈ (2.5, 5.0)) + @test all(inv_spacing(grid, (Vertex(), Center()), 1, 1) .≈ (2.5, 5.0)) + # repeating locations + @test inv_spacing(grid, Vertex(), 1, 1) == inv_spacing(grid, (Vertex(), Vertex()), 1, 1) + @test inv_spacing(grid, Center(), 1, 1) == inv_spacing(grid, (Center(), Center()), 1, 1) + end + + @testset "coords" begin + # one index + @test coord(grid, Vertex(), Dim(1), 1) ≈ -1.0 + @test coord(grid, Vertex(), Dim(2), 1) ≈ -2.0 + @test coord(grid, Vertex(), Dim(1), nx + 1) ≈ 1.0 + @test coord(grid, Vertex(), Dim(2), ny + 1) ≈ 2.0 + @test coord(grid, Center(), Dim(1), 1) ≈ -0.8 + @test coord(grid, Center(), Dim(2), 1) ≈ -1.9 + @test coord(grid, Center(), Dim(1), nx) ≈ 0.8 + @test coord(grid, Center(), Dim(2), ny) ≈ 1.9 + # many indices + for loc in (Center(), Vertex()) + @test coord(grid, loc, Dim(1), 1, 1) == coord(grid, loc, Dim(1), 1) + @test coord(grid, loc, Dim(2), 1, 1) == coord(grid, loc, Dim(2), 1) + @test coord(grid, loc, Dim(1), nx + 1, 1) == coord(grid, loc, Dim(1), nx + 1) + @test coord(grid, loc, Dim(2), 1, ny + 1) == coord(grid, loc, Dim(2), ny + 1) + end + # many locations + @test all(coord(grid, (Vertex(), Center()), 1, 1) .≈ (-1.0, -1.9)) + @test all(coord(grid, (Vertex(), Center()), nx + 1, 1) .≈ (1.0, -1.9)) + @test all(coord(grid, (Vertex(), Center()), 1, ny) .≈ (-1.0, 1.9)) + # repeated locations + @test coord(grid, (Vertex(), Center()), Dim(1), 1) == coord(grid, Vertex(), Dim(1), 1) == coord(grid, Vertex(), Dim(1), 1, 1) + @test coord(grid, (Vertex(), Center()), Dim(2), 1) == coord(grid, Center(), Dim(2), 1) == coord(grid, Center(), Dim(2), 1, 1) + end + + @testset "shortcut coords" begin + # short coords + @test vertex(grid, Dim(1), 1) == vertex(grid, Dim(1), 1, 1) == coord(grid, Vertex(), Dim(1), 1) + @test vertex(grid, Dim(2), 1) == vertex(grid, Dim(2), 1, 1) == coord(grid, Vertex(), Dim(2), 1) + @test center(grid, Dim(1), 1) == center(grid, Dim(1), 1, 1) == coord(grid, Center(), Dim(1), 1) + @test center(grid, Dim(2), 1) == center(grid, Dim(2), 1, 1) == coord(grid, Center(), Dim(2), 1) + end + + @testset "cartesian" begin + # coords + @test xvertex(grid, 1) == vertex(grid, Dim(1), 1) + @test xcenter(grid, 1) == center(grid, Dim(1), 1) + @test yvertex(grid, 1) == vertex(grid, Dim(2), 1) + @test ycenter(grid, 1) == center(grid, Dim(2), 1) + # spacing + for loc in (Center(), Vertex()) + @test Δx(grid, loc, 1) == spacing(grid, loc, Dim(1), 1) + @test Δx(grid, loc, 1) == spacing(grid, loc, Dim(1), 1) + @test Δy(grid, loc, 1) == spacing(grid, loc, Dim(2), 1) + @test Δy(grid, loc, 1) == spacing(grid, loc, Dim(2), 1) + end + end + end + end +end diff --git a/test/test_distributed_2D.jl b/test/test_distributed_2D.jl deleted file mode 100644 index c63084f2..00000000 --- a/test/test_distributed_2D.jl +++ /dev/null @@ -1,49 +0,0 @@ -include("common.jl") - -using MPI -using FastIce.Distributed - -MPI.Init() - -nprocs = MPI.Comm_size(MPI.COMM_WORLD) - -backends = [CPU()] # until we have testing environment setup for GPU-aware MPI, run only on CPU - -for backend in backends - @testset "$(basename(@__FILE__)) (backend: $backend)" begin - mpi_dims = parse.(Int, (get(ENV, "DIMX", "1"), - get(ENV, "DIMY", "1"))) - dims = (0, 0) - topo = CartesianTopology(dims) - local_size = (4, 5) - @testset "Topology" begin - @test dimensions(topo) == mpi_dims - @test length(neighbors(topo)) == 2 - @test node_size(topo) == nprocs - if global_rank(topo) == 0 - @test neighbor(topo, 1, 1) == MPI.PROC_NULL - @test neighbor(topo, 2, 1) == MPI.PROC_NULL - @test has_neighbor(topo, 1, 1) == false - @test has_neighbor(topo, 2, 1) == false - if mpi_dims[2] > 1 - @test neighbor(topo, 1, 2) == mpi_dims[2] - @test has_neighbor(topo, 1, 2) == true - else - @test neighbor(topo, 1, 2) == MPI.PROC_NULL - @test has_neighbor(topo, 1, 2) == false - end - end - @test global_grid_size(topo, local_size) == mpi_dims .* local_size - end - @testset "gather!" begin - src = fill(global_rank(topo) + 1, local_size) - dst = (global_rank(topo) == 0) ? zeros(Int, mpi_dims .* local_size) : nothing - gather!(dst, src, cartesian_communicator(topo)) - if global_rank(topo) == 0 - @test dst == repeat(reshape(1:global_size(topo), dimensions(topo))'; inner=size(src)) - end - end - end -end - -MPI.Finalize() diff --git a/test/test_distributed_3D.jl b/test/test_distributed_3D.jl deleted file mode 100644 index a5a24841..00000000 --- a/test/test_distributed_3D.jl +++ /dev/null @@ -1,57 +0,0 @@ -include("common.jl") - -using MPI -using FastIce.Distributed - -MPI.Init() - -nprocs = MPI.Comm_size(MPI.COMM_WORLD) - -backends = [CPU()] # until we have testing environment setup for GPU-aware MPI, run only on CPU - -for backend in backends - @testset "$(basename(@__FILE__)) (backend: $backend)" begin - mpi_dims = parse.(Int, (get(ENV, "DIMX", "1"), - get(ENV, "DIMY", "1"), - get(ENV, "DIMZ", "1"))) - dims = (0, 0, 0) - topo = CartesianTopology(dims) - local_size = (4, 5, 6) - @testset "Topology" begin - @test dimensions(topo) == mpi_dims - @test length(neighbors(topo)) == 3 - @test node_size(topo) == nprocs - if global_rank(topo) == 0 - @test neighbor(topo, 1, 1) == MPI.PROC_NULL - @test neighbor(topo, 2, 1) == MPI.PROC_NULL - @test neighbor(topo, 3, 1) == MPI.PROC_NULL - @test has_neighbor(topo, 1, 1) == false - @test has_neighbor(topo, 2, 1) == false - @test has_neighbor(topo, 3, 1) == false - if mpi_dims[2] > 1 && mpi_dims[3] > 1 - @test neighbor(topo, 1, 2) == mpi_dims[2] * mpi_dims[3] - @test neighbor(topo, 2, 2) == mpi_dims[3] - @test has_neighbor(topo, 1, 2) == true - @test has_neighbor(topo, 2, 2) == true - else - @test neighbor(topo, 1, 2) == MPI.PROC_NULL - @test neighbor(topo, 2, 2) == MPI.PROC_NULL - @test has_neighbor(topo, 1, 2) == false - @test has_neighbor(topo, 2, 2) == false - end - end - @test global_grid_size(topo, local_size) == mpi_dims .* local_size - end - @testset "gather!" begin - src = fill(global_rank(topo) + 1, local_size) - dst = (global_rank(topo) == 0) ? zeros(Int, mpi_dims .* local_size) : nothing - gather!(dst, src, cartesian_communicator(topo)) - ranks_mat = permutedims(reshape(1:global_size(topo), dimensions(topo)), reverse(1:3)) - if global_rank(topo) == 0 - @test dst == repeat(ranks_mat; inner=size(src)) - end - end - end -end - -MPI.Finalize() diff --git a/test/test_distributed_writers_3D.jl b/test/test_distributed_writers_3D.jl index 22868b02..9e6e5314 100644 --- a/test/test_distributed_writers_3D.jl +++ b/test/test_distributed_writers_3D.jl @@ -1,9 +1,10 @@ include("common.jl") -using FastIce.Architectures -using FastIce.Distributed -using FastIce.Fields -using FastIce.Grids +using Chmy.Architectures +using Chmy.Distributed +using Chmy.Fields +using Chmy.Grids + using FastIce.Writers using MPI @@ -14,13 +15,7 @@ MPI.Init() backends = [CPU()] # until we have testing environment setup for GPU-aware MPI, run only on CPU -dims = (0, 0, 0) -topo = CartesianTopology(dims) -mpi_dims = dimensions(topo) -me = global_rank(topo) # rank -comm = cartesian_communicator(topo) - -(me == 0) && (XML_ref = """ +XML_ref = """ @@ -29,43 +24,47 @@ comm = cartesian_communicator(topo) -""") +""" +dims_g = (0, 0, 0) size_l = (4, 5, 6) -size_g = global_grid_size(topo, size_l) -Fa_g = (me == 0) ? zeros(Float64, mpi_dims .* size_l) : nothing -Fb_g = (me == 0) ? zeros(Float64, mpi_dims .* size_l) : nothing +h5_fname = "test_d.h5" +xdmf_fname = "test_d.xdmf3" -grid_g = CartesianGrid(; origin=(-0.4, -0.5, 0.0), - extent=(1.0, 1.1, 1.2), - size=size_g) +for backend in backends -grid_l = local_grid(grid_g, topo) + arch = Arch(backend, MPI.COMM_WORLD, dims_g) + topo = topology(arch) + me = global_rank(topo) # rank + mpi_dims = dims(topo) + comm = cart_comm(topo) + + size_g = size_l .* dims(topo) + grid = UniformGrid(arch; origin=(-0.4, -0.5, 0.0), extent=(1.0, 1.1, 1.2), dims=size_g) + + Fa_g = (me == 0) ? zeros(Float64, size_g) : nothing + Fb_g = (me == 0) ? zeros(Float64, size_g) : nothing -for backend in backends @testset "$(basename(@__FILE__)) (backend: $backend)" begin HDF5.has_parallel() || (@warn("HDF5 has no parallel support. Skipping $(basename(@__FILE__)) (backend: $backend)."); return) - arch = Architecture(backend, topo) - set_device!(arch) - - Fa_l = Field(backend, grid_l, Center()) - Fb_l = Field(backend, grid_l, Center()) + Fa_l = Field(backend, grid, Center()) + Fb_l = Field(backend, grid, Center()) fill!(parent(Fa_l), 1.0 * me) fill!(parent(Fb_l), 2.0 * me) @@ -73,37 +72,30 @@ for backend in backends fields = Dict("Fa" => Fa_l, "Fb" => Fb_l) @testset "Distributed writers 3D" begin - @testset "write/read h5" begin - fname = "test_d.h5" - (me == 0) && (isfile(fname) && run(`rm $fname`)) - write_h5(arch, grid_g, fname, fields) + @testset "write h5" begin + (me == 0) && (isfile(h5_fname) && run(`rm $h5_fname`)) + write_h5(arch, grid, h5_fname, fields) MPI.Barrier(comm) - Fa_v = zeros(size(grid_l)) - Fb_v = zeros(size(grid_l)) - copyto!(Fa_v, interior(Fa_l)) - copyto!(Fb_v, interior(Fb_l)) - KernelAbstractions.synchronize(backend) - - gather!(Fa_g, Fa_v, comm) - gather!(Fb_g, Fb_v, comm) + gather!(arch, Fa_g, Fa_l) + gather!(arch, Fb_g, Fb_l) MPI.Barrier(comm) if me == 0 - @test all(Fa_g .== h5read(fname, "Fa")) - @test all(Fb_g .== h5read(fname, "Fb")) - isfile(fname) && run(`rm $fname`) + @test all(Fa_g .== h5read(h5_fname, "Fa")) + @test all(Fb_g .== h5read(h5_fname, "Fb")) + isfile(h5_fname) && run(`rm $h5_fname`) end end MPI.Barrier(comm) + @testset "write/read xdmf3" begin if me == 0 - fname = "test_d.xdmf3" - isfile(fname) && run(`rm $fname`) - write_xdmf(arch, grid_g, fname, fields, "test_d.h5") + isfile(xdmf_fname) && run(`rm $xdmf_fname`) + write_xdmf(arch, grid, xdmf_fname, fields, (h5_fname, )) - @test XML_ref == string(parse_file(fname)) - isfile(fname) && run(`rm $fname`) + @test XML_ref == string(parse_file(xdmf_fname)) + isfile(xdmf_fname) && run(`rm $xdmf_fname`) end end end diff --git a/test/test_fields.jl b/test/test_fields.jl deleted file mode 100644 index 0e93de1a..00000000 --- a/test/test_fields.jl +++ /dev/null @@ -1,61 +0,0 @@ -include("common.jl") - -using FastIce.Fields -using FastIce.Grids - -const HALOS = ((1, 1, 2), nothing, (1, nothing, 1), (nothing, nothing, nothing), (nothing, (1, nothing), (1, 2))) -const HALO_SIZES = ((2, 2, 4), (0, 0, 0), (2, 0, 2), (0, 0, 0), (0, 1, 3)) - -for backend in backends - @testset "$(basename(@__FILE__)) (backend: $backend)" begin - grid = CartesianGrid(; origin=(0.0, 0.0, 0.0), extent=(1.0, 1.0, 1.0), size=(2, 2, 2)) - loc = (Center(), Vertex(), Center()) - @testset "location" begin - @test location(Field(backend, grid, Center())) == (Center(), Center(), Center()) - @test location(Field(backend, grid, loc)) == loc - end - @testset "halo $hl" for (hl, hs) in zip(HALOS, HALO_SIZES) - f = Field(backend, grid, loc; halo=hl) - @test location(f) == (Center(), Vertex(), Center()) - @test size(data(f)) == size(grid, loc) .+ hs - @test size(interior(f)) == size(grid, loc) - end - @testset "set" begin - f = Field(backend, grid, (Center(), Vertex(), Center()); halo=(1, 0, 1)) - @testset "discrete" begin - # no parameters vertex - fill!(data(f), NaN) - set!(f, grid, (grid, loc, I) -> ycoord(grid, loc, I[2]); discrete=true) - @test Array(interior(f)) == [0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0;;; - 0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0] - # no parameters center - fill!(data(f), NaN) - set!(f, grid, (grid, loc, I) -> xcoord(grid, loc, I[1]); discrete=true) - @test Array(interior(f)) == [0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75;;; - 0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75] - # with parameters - fill!(data(f), NaN) - set!(f, grid, (grid, loc, I, sc) -> ycoord(grid, loc, I[2]) * sc; discrete=true, parameters=(2.0,)) - @test Array(interior(f)) == [0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0;;; - 0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0] - end - @testset "continuous" begin - # no parameters vertex - fill!(data(f), NaN) - set!(f, grid, (x, y, z) -> y) - @test Array(interior(f)) == [0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0;;; - 0.0; 0.0;; 0.5; 0.5;; 1.0; 1.0] - # no parameters center - fill!(data(f), NaN) - set!(f, grid, (x, y, z) -> x) - @test Array(interior(f)) == [0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75;;; - 0.25; 0.75;; 0.25; 0.75;; 0.25; 0.75] - # with parameters - fill!(data(f), NaN) - set!(f, grid, (x, y, z, sc) -> y * sc; parameters=(2.0,)) - @test Array(interior(f)) == [0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0;;; - 0.0; 0.0;; 1.0; 1.0;; 2.0; 2.0] - end - end - end -end diff --git a/test/test_utils.jl b/test/test_utils.jl deleted file mode 100644 index 29c9d767..00000000 --- a/test/test_utils.jl +++ /dev/null @@ -1,41 +0,0 @@ -include("common.jl") - -using FastIce.Utils - -@testset "$(basename(@__FILE__)) (backend: CPU)" begin - @testset "workers" begin - @testset "setup" begin - a = 0 - worker = Worker(; setup=() -> a += 1) - put!(() -> nothing, worker) - wait(worker) - @test a == 1 - close(worker) - end - @testset "teardown" begin - a = 0 - worker = Worker(; teardown=() -> a += 2) - put!(worker) do - a -= 1 - end - wait(worker) - close(worker) - @test a == 1 - end - @testset "do work" begin - a = 0 - worker = Worker() - put!(worker) do - a += 1 - end - wait(worker) - close(worker) - @test a == 1 - end - @testset "not running" begin - worker = Worker() - close(worker) - @test_throws ErrorException wait(worker) - end - end -end diff --git a/test/test_writers.jl b/test/test_writers.jl index a86908f7..904ee3a0 100644 --- a/test/test_writers.jl +++ b/test/test_writers.jl @@ -1,8 +1,9 @@ include("common.jl") -using FastIce.Architectures -using FastIce.Fields -using FastIce.Grids +using Chmy.Architectures +using Chmy.Fields +using Chmy.Grids + using FastIce.Writers using HDF5 @@ -17,14 +18,14 @@ XML_ref = """