From 58eb20e91c69b9eef973fdf2394dc5382ce63e07 Mon Sep 17 00:00:00 2001 From: Ivan Utkin Date: Sun, 3 Mar 2024 13:24:28 +0100 Subject: [PATCH 01/71] Port to Chmy (WIP) --- Project.toml | 11 +- ext/FastIceAMDGPUExt/FastIceAMDGPUExt.jl | 14 -- ext/FastIceCUDAExt/FastIceCUDAExt.jl | 15 -- scripts_future_API/bench3d.jl | 139 ------------- ...okes_manufactured_solution_free_slip_2D.jl | 32 +-- ...okes_manufactured_solution_free_slip_3D.jl | 17 +- src/Architectures.jl | 43 ---- src/BoundaryConditions/BoundaryConditions.jl | 63 ------ src/BoundaryConditions/boundary_function.jl | 39 ---- src/BoundaryConditions/dirichlet_bc.jl | 34 --- src/BoundaryConditions/extrapolate_bc.jl | 29 --- .../field_boundary_conditions.jl | 44 ---- src/BoundaryConditions/hide_boundaries.jl | 30 --- src/BoundaryConditions/utils.jl | 23 --- src/Distributed/Distributed.jl | 38 ---- src/Distributed/boundary_conditions.jl | 112 ---------- src/Distributed/gather.jl | 20 -- src/Distributed/topology.jl | 157 -------------- src/FastIce.jl | 8 - src/Fields/Fields.jl | 156 -------------- src/Fields/constant_field.jl | 16 -- src/Fields/function_field.jl | 46 ----- src/GridOperators.jl | 77 ------- src/Grids/Grids.jl | 18 -- src/Grids/cartesian_grid.jl | 193 ------------------ src/Grids/discrete_axis.jl | 61 ------ src/KernelLaunch.jl | 86 -------- src/LevelSets/LevelSets.jl | 10 +- src/LevelSets/compute_level_sets.jl | 14 +- src/LevelSets/signed_distances.jl | 4 +- src/Models/Diffusion/Heat/Heat.jl | 71 +++---- .../Diffusion/Heat/boundary_conditions.jl | 14 -- .../FullStokes/Isothermal/Isothermal.jl | 126 ++++-------- .../FullStokes/Isothermal/kernels_2d.jl | 99 ++++----- .../FullStokes/Isothermal/kernels_3d.jl | 152 ++++++-------- .../FullStokes/Isothermal/kernels_common.jl | 15 -- src/Utils/Utils.jl | 36 ---- src/Utils/split_ndrange.jl | 26 --- src/Utils/workers.jl | 45 ---- src/Writers.jl | 57 +++--- test/test_bc.jl | 93 --------- test/test_distributed_2D.jl | 49 ----- test/test_distributed_3D.jl | 57 ------ test/test_fields.jl | 61 ------ test/test_utils.jl | 41 ---- test/test_writers.jl | 8 +- 46 files changed, 227 insertions(+), 2272 deletions(-) delete mode 100644 ext/FastIceAMDGPUExt/FastIceAMDGPUExt.jl delete mode 100644 ext/FastIceCUDAExt/FastIceCUDAExt.jl delete mode 100644 scripts_future_API/bench3d.jl delete mode 100644 src/Architectures.jl delete mode 100644 src/BoundaryConditions/BoundaryConditions.jl delete mode 100644 src/BoundaryConditions/boundary_function.jl delete mode 100644 src/BoundaryConditions/dirichlet_bc.jl delete mode 100644 src/BoundaryConditions/extrapolate_bc.jl delete mode 100644 src/BoundaryConditions/field_boundary_conditions.jl delete mode 100644 src/BoundaryConditions/hide_boundaries.jl delete mode 100644 src/BoundaryConditions/utils.jl delete mode 100644 src/Distributed/Distributed.jl delete mode 100644 src/Distributed/boundary_conditions.jl delete mode 100644 src/Distributed/gather.jl delete mode 100644 src/Distributed/topology.jl delete mode 100644 src/Fields/Fields.jl delete mode 100644 src/Fields/constant_field.jl delete mode 100644 src/Fields/function_field.jl delete mode 100644 src/GridOperators.jl delete mode 100644 src/Grids/Grids.jl delete mode 100644 src/Grids/cartesian_grid.jl delete mode 100644 src/Grids/discrete_axis.jl delete mode 100644 src/KernelLaunch.jl delete mode 100644 src/Models/FullStokes/Isothermal/kernels_common.jl delete mode 100644 src/Utils/Utils.jl delete mode 100644 src/Utils/split_ndrange.jl delete mode 100644 src/Utils/workers.jl delete mode 100644 test/test_bc.jl delete mode 100644 test/test_distributed_2D.jl delete mode 100644 test/test_distributed_3D.jl delete mode 100644 test/test_fields.jl delete mode 100644 test/test_utils.jl diff --git a/Project.toml b/Project.toml index 6d5a66d8..7472a4f5 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.2.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" @@ -18,18 +19,8 @@ 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" 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/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/stokes_manufactured_solution_free_slip_2D.jl b/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl index 07e24c4d..5efbf3e1 100644 --- a/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl +++ b/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl @@ -1,14 +1,14 @@ using Printf -using FastIce -using FastIce.Architectures -using FastIce.Grids -using FastIce.Fields -using FastIce.Utils -using FastIce.BoundaryConditions +using Chmy +using Chmy.Architectures +using Chmy.Grids +using Chmy.Fields +using Chmy.BoundaryConditions +using Chmy.Physics +using Chmy.KernelLaunch + using FastIce.Models.FullStokes.Isothermal -using FastIce.Physics -using FastIce.KernelLaunch const VBC = BoundaryCondition{Velocity} const TBC = BoundaryCondition{Traction} @@ -38,7 +38,7 @@ vy(x, y) = -cos(π * x) * sin(π * y) @views function run(dims) backend = CPU() - arch = Architecture(backend, 2) + arch = Arch(backend) set_device!(arch) # physics @@ -46,9 +46,9 @@ vy(x, y) = -cos(π * x) * sin(π * y) A0 = 0.5 # geometry - grid = CartesianGrid(; origin=(-1.0, -1.0), - extent=(2.0, 2.0), - size=dims) + grid = UniformGrid(; origin=(-1.0, -1.0), + extent=(2.0, 2.0), + dims) free_slip = SBC(0.0, 0.0) xface = (Vertex(), Center()) @@ -93,10 +93,10 @@ vy(x, y) = -cos(π * x) * sin(π * y) set!(model.viscosity.η, 1.0) set!(model.viscosity.η_next, 1.0) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.stress) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.velocity) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.viscosity.η) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.viscosity.η_next) + bc!(arch, grid, model.boundary_conditions.stress) + bc!(arch, grid, model.boundary_conditions.velocity) + bc!(arch, grid, model.boundary_conditions.viscosity.η) + bc!(arch, grid, model.boundary_conditions.viscosity.η_next) for iter in 1:niter advance_iteration!(model, 0.0, 1.0) diff --git a/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl b/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl index 1f09ec77..5d43daec 100644 --- a/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl +++ b/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl @@ -1,14 +1,15 @@ using Printf -using FastIce -using FastIce.Architectures -using FastIce.Grids -using FastIce.Fields -using FastIce.Utils -using FastIce.BoundaryConditions +using Chmy +using Chmy.Architectures +using Chmy.Grids +using Chmy.Fields +using Chmy.Utils +using Chmy.BoundaryConditions +using Chmy.Physics +using Chmy.KernelLaunch + using FastIce.Models.FullStokes.Isothermal -using FastIce.Physics -using FastIce.KernelLaunch const VBC = BoundaryCondition{Velocity} const TBC = BoundaryCondition{Traction} 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..143c7725 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -26,15 +26,7 @@ greet(; kwargs...) = printstyled(GREETING; kwargs...) greet_fast(; kwargs...) = printstyled(GREETING_FAST; kwargs...) # core modules -include("Grids/Grids.jl") -include("GridOperators.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") 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/GridOperators.jl b/src/GridOperators.jl deleted file mode 100644 index 5b542a69..00000000 --- a/src/GridOperators.jl +++ /dev/null @@ -1,77 +0,0 @@ -module GridOperators - -export ∂ᶜ, ∂ᵛ - -export ∂ᵛx, ∂ᵛy, ∂ᵛz, avᵛx, avᵛy, avᵛz, avᵛxy, avᵛxz, avᵛyz, maxlᵛx, maxlᵛy, maxlᵛz -export ∂ᶜx, ∂ᶜy, ∂ᶜz, avᶜx, avᶜy, avᶜz, avᶜxy, avᶜxz, avᶜyz, maxlᶜx, maxlᶜy, maxlᶜz - -export maxlᵛxy, maxlᵛxz, maxlᵛyz -export maxlᵛyx, maxlᵛzx, maxlᵛzy - -import Base.@propagate_inbounds - -Base.@assume_effects :foldable function δ(op, I::CartesianIndex{N}, ::Val{D}) where {N,D} - ntuple(i -> 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..5eb8428c 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -2,14 +2,14 @@ module LevelSets export compute_level_set_from_dem! -using FastIce -using FastIce.Grids -using FastIce.Architectures -using FastIce.Fields +using Chmy +using Chmy.Grids +using Chmy.Architectures +using Chmy.Fields using KernelAbstractions using LinearAlgebra, GeometryBasics include("signed_distances.jl") include("compute_level_sets.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 index ccd2d26b..3bade50b 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -1,9 +1,4 @@ -""" - _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) +@kernel inbounds = true function _init_level_set!(Ψ::Field, dem::Field, dem_grid::StructuredGrid, Ψ_grid::StructuredGrid, cutoff, R) I = @index(Global, Cartesian) x, y, z = coord(Ψ_grid, location(Ψ), I) P = R * Point3(x, y, z) @@ -12,11 +7,12 @@ Initialize level sets. end """ - compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, R=LinearAlgebra.I) + compute_level_set_from_dem!(arch, Ψ, dem, dem_grid, Ψ_grid, R[=LinearAlgebra.I]) -Compute level sets from dem. +Compute a level set from a given dem. """ -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, R=LinearAlgebra.I) +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) diff --git a/src/LevelSets/signed_distances.jl b/src/LevelSets/signed_distances.jl index 5e4b9da6..7d4c43f9 100644 --- a/src/LevelSets/signed_distances.jl +++ b/src/LevelSets/signed_distances.jl @@ -1,5 +1,3 @@ -export sd_dem - @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))) @@ -69,4 +67,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/Models/Diffusion/Heat/Heat.jl b/src/Models/Diffusion/Heat/Heat.jl index 967bd79c..8f97bb8c 100644 --- a/src/Models/Diffusion/Heat/Heat.jl +++ b/src/Models/Diffusion/Heat/Heat.jl @@ -3,75 +3,56 @@ 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 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/FullStokes/Isothermal/Isothermal.jl b/src/Models/FullStokes/Isothermal/Isothermal.jl index 5c951a61..b6ee2ede 100644 --- a/src/Models/FullStokes/Isothermal/Isothermal.jl +++ b/src/Models/FullStokes/Isothermal/Isothermal.jl @@ -3,16 +3,16 @@ 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.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") @@ -21,7 +21,7 @@ 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,56 +31,25 @@ 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())))) +function StressFields(arch, grid) + return (Pr=Field(arch, grid, Center()), τ=TensorField(arch, grid)) end function VelocityFields(arch, grid::CartesianGrid{2}) - (x=Field(arch, grid, (Vertex(), Center()); halo=1), - y=Field(arch, grid, (Center(), Vertex()); halo=1)) + return VectorField(arch, grid) 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())))) -end - -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())))) -end - -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)) -end - -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())))) + (r_Pr=Field(arch, grid, Center()), r_V=VectorField(arch, grid)) end function ViscosityFields(arch, grid::CartesianGrid) - (η = Field(arch, grid, Center(); halo=1), - η_next = Field(arch, grid, Center(); halo=1)) + (η = Field(arch, grid, Center()), + η_next = Field(arch, grid, Center())) end function IsothermalFullStokesModel(; arch, @@ -101,7 +70,7 @@ function IsothermalFullStokesModel(; arch, outer_width = ntuple(_ -> 2, Val(ndims(grid))) end - hide_boundaries = HideBoundaries(arch, outer_width) + launcher = Launcher(arch, grid; outer_width) return IsothermalFullStokesModel(arch, grid, @@ -113,56 +82,39 @@ function IsothermalFullStokesModel(; arch, 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) + (; Pr, τ) = model.stress + V = model.velocity + (; η, η_next) = model.viscosity + (; Δτ) = model.solver_params + arch = model.arch + rheology = model.rheology + ρg = model.gravity + bc = model.boundary_conditions + grid = model.grid + launch = model.launcher + + launch(arch, grid, update_σ! => (Pr, τ, V, η, Δτ, grid); bc=bc.stress) + launch(arch, grid, update_V! => (V, Pr, τ, η, η_next, rheology, ρg, Δτ, grid); bc=velocity_bc) # 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 + (; 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 + grid = model.grid + ρg = model.gravity + bc = model.boundary_conditions.residual + launch!(model.arch, grid, compute_residuals! => (r_V, r_Pr, Pr, τ, V, ρg, grid); bc) return end diff --git a/src/Models/FullStokes/Isothermal/kernels_2d.jl b/src/Models/FullStokes/Isothermal/kernels_2d.jl index f13a8e24..55230281 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 * η[I] * Δτ.Pr + # deviatoric stress + τ.xx[I] -= (τ.xx[I] - 2.0 * itp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx + τ.yy[I] -= (τ.yy[I] - 2.0 * itp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy + τ.xy[I] -= (τ.xy[I] - 2.0 * itp(η, 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]) / itp(η, location(V.x), I) * Δτ.V.x + V.y[I] += (-∂y(Pr, g, I) + ∂y(τ.yy, g, I) + ∂x(τ.xy, g, I) - ρg.y[I]) / itp(η, location(V.x), I) * Δτ.V.y + # rheology + τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2) + itp(τ.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::CartesianGrid{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 * itp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0) + τ.yy[I] = 2.0 * itp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0) + τ.xy[I] = 2.0 * itp(η, 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::CartesianGrid{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..98f8044e 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::CartesianGrid{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 * η[I] * Δτ.Pr + # deviatoric diagonal + τ.xx[I] -= (τ.xx[I] - 2.0 * itp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0)) * Δτ.τ.xx + τ.yy[I] -= (τ.yy[I] - 2.0 * itp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0)) * Δτ.τ.yy + τ.zz[I] -= (τ.zz[I] - 2.0 * itp(η, 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 * itp(η, location(τ.xy), g, I) * ε̇xy) * Δτ.τ.xy + τ.xz[I] -= (τ.xz[I] - 2.0 * itp(η, location(τ.xz), g, I) * ε̇xz) * Δτ.τ.xz + τ.yz[I] -= (τ.yz[I] - 2.0 * itp(η, 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::CartesianGrid{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]) / itp(η, 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]) / itp(η, 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]) / itp(η, location(V.z), g, I) * Δτ.V.z + # rheology + τII = sqrt(0.5 * (τ.xx[I]^2 + τ.yy[I]^2 + τ.zz[I]^2) + + itp(τ.xy, location(η), g, I)^2 + + itp(τ.xz, location(η), g, I)^2 + + itp(τ.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::CartesianGrid{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 * itp(η, location(τ.xx), g, I) * (ε̇xx - ∇V / 3.0) + τ.yy[I] = 2.0 * itp(η, location(τ.yy), g, I) * (ε̇yy - ∇V / 3.0) + τ.zz[I] = 2.0 * itp(η, location(τ.zz), g, I) * (ε̇zz - ∇V / 3.0) + # deviatoric off-diagonal + τ.xy[I] = 2.0 * itp(η, location(τ.xy), g, I) * ε̇xy + τ.xz[I] = 2.0 * itp(η, location(τ.xz), g, I) * ε̇xz + τ.yz[I] = 2.0 * itp(η, 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::CartesianGrid{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/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..f97204ac 100644 --- a/src/Writers.jl +++ b/src/Writers.jl @@ -2,23 +2,31 @@ 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::CartesianGrid, path, fields) + I = CartesianIndices(size(grid)) + h5open(path, "w") do io + write_dset(io, fields, size(grid), I.indices) + end + return +end + +function write_h5(arch::DistributedArchitecture, grid::CartesianGrid, path, fields) HDF5.has_parallel() || @warn("HDF5 has no parallel support.") - topo = details(arch) + topo = topology(arch) comm = cartesian_communicator(topo) coords = coordinates(topo) sz = size(local_grid(grid, topo)) @@ -31,19 +39,6 @@ function write_h5(arch::Architecture{DistributedMPI}, grid::CartesianGrid, path, 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) - end - return -end - function write_dset(io, fields, grid_size, inds) for (name, field) in fields dset = create_dataset(io, "/$name", datatype(eltype(field)), dataspace(grid_size)) @@ -53,16 +48,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) +function write_xdmf(arch::Architecture, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) grid_size = size(grid) grid_spacing = spacing(grid) - grid_origin = origin(local_grid(grid, topo)) + grid_origin = origin(grid) xdoc = generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) @@ -70,16 +64,11 @@ 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)) +function write_xdmf(arch::DistributedArchitecture, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) + topo = topology(arch) grid_size = size(grid) grid_spacing = spacing(grid) - grid_origin = origin(grid) + grid_origin = origin(local_grid(grid, topo)) xdoc = generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) 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_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_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..f64d9b50 100644 --- a/test/test_writers.jl +++ b/test/test_writers.jl @@ -1,8 +1,8 @@ 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 @@ -38,7 +38,7 @@ grid = CartesianGrid(; origin=(-0.4, -0.5, 0.0), for backend in backends @testset "$(basename(@__FILE__)) (backend: $backend)" begin - arch = Architecture(backend) + arch = Arch(backend) Fa = Field(backend, grid, Center()) Fb = Field(backend, grid, Center()) From 6f160714f08f0c254d7ee17f096cf74a38bbc289 Mon Sep 17 00:00:00 2001 From: Ivan Utkin Date: Sun, 3 Mar 2024 13:28:23 +0100 Subject: [PATCH 02/71] Port to Chmy (WIP) --- ...okes_manufactured_solution_free_slip_2D.jl | 8 +++--- ...okes_manufactured_solution_free_slip_3D.jl | 24 +++++++++--------- ...s_manufactured_solution_free_surface_2D.jl | 25 ++++++++++--------- 3 files changed, 29 insertions(+), 28 deletions(-) diff --git a/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl b/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl index 5efbf3e1..9e7a58c8 100644 --- a/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl +++ b/scripts_future_API/stokes_manufactured_solution_free_slip_2D.jl @@ -61,14 +61,14 @@ vy(x, y) = -cos(π * x) * sin(π * y) y=FunctionField(ρgy, grid, yface; parameters=η0)) # numerics - niter = 10maximum(size(grid)) - ncheck = maximum(size(grid)) + niter = 10maximum(size(grid, Center())) + ncheck = maximum(size(grid, Center())) # PT params r = 0.7 re_mech = 8π - lτ_re_m = minimum(extent(grid)) / re_mech - vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.1) + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid, Center(), 1, 1)) / 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 diff --git a/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl b/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl index 5d43daec..e3a2c7f4 100644 --- a/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl +++ b/scripts_future_API/stokes_manufactured_solution_free_slip_3D.jl @@ -57,7 +57,7 @@ vz(x, y, z) = sin(π * z) * f(x, y) @views function run(dims) backend = CUDABackend() - arch = Architecture(backend, 1) + arch = Arch(backend) set_device!(arch) # outer_width = (4, 4, 4) @@ -67,9 +67,9 @@ vz(x, y, z) = sin(π * z) * f(x, y) A0 = 0.5 # geometry - grid = CartesianGrid(; origin=(-1.0, -1.0, -1.0), - extent=(2.0, 2.0, 2.0), - size=dims) + grid = UniformGrid(; origin=(-1.0, -1.0, -1.0), + extent=(2.0, 2.0, 2.0), + dims) free_slip = SBC(0.0, 0.0, 0.0) xface = (Vertex(), Center(), Center()) @@ -85,14 +85,14 @@ vz(x, y, z) = sin(π * z) * f(x, y) z=FunctionField(ρgz, grid, zface; parameters=η0)) # numerics - niter = 10maximum(size(grid)) - ncheck = maximum(size(grid)) + niter = 10maximum(size(grid, Center())) + ncheck = maximum(size(grid, Center())) # PT params r = 0.7 re_mech = 8π - lτ_re_m = minimum(extent(grid)) / re_mech - vdτ = minimum(spacing(grid)) / sqrt(ndims(grid) * 1.1) + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid, Center(), 1, 1)) / 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 @@ -117,10 +117,10 @@ vz(x, y, z) = sin(π * z) * f(x, y) set!(model.viscosity.η, 1.0) set!(model.viscosity.η_next, 1.0) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.stress) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.velocity) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.viscosity.η) - KernelLaunch.apply_all_boundary_conditions!(arch, grid, model.boundary_conditions.viscosity.η_next) + bc!(arch, grid, model.boundary_conditions.stress) + bc!(arch, grid, model.boundary_conditions.velocity) + bc!(arch, grid, model.boundary_conditions.viscosity.η) + bc!(arch, grid, model.boundary_conditions.viscosity.η_next) for iter in 1:niter advance_iteration!(model, 0.0, 1.0) diff --git a/scripts_future_API/stokes_manufactured_solution_free_surface_2D.jl b/scripts_future_API/stokes_manufactured_solution_free_surface_2D.jl index a7820bd0..3b5cb6a3 100644 --- a/scripts_future_API/stokes_manufactured_solution_free_surface_2D.jl +++ b/scripts_future_API/stokes_manufactured_solution_free_surface_2D.jl @@ -1,14 +1,15 @@ using Printf -using FastIce -using FastIce.Architectures -using FastIce.Grids -using FastIce.Fields -using FastIce.Utils -using FastIce.BoundaryConditions +using Chmy +using Chmy.Architectures +using Chmy.Grids +using Chmy.Fields +using Chmy.Utils +using Chmy.BoundaryConditions +using Chmy.Physics +using Chmy.KernelLaunch + using FastIce.Models.FullStokes.Isothermal -using FastIce.Physics -using FastIce.KernelLaunch const VBC = BoundaryCondition{Velocity} const TBC = BoundaryCondition{Traction} @@ -38,7 +39,7 @@ vy(x, y) = -cos(0.5π * x) * sin(0.5π * y) @views function run(dims) backend = CPU() - arch = Architecture(backend, 2) + arch = Arch(backend) set_device!(arch) # physics @@ -48,7 +49,7 @@ vy(x, y) = -cos(0.5π * x) * sin(0.5π * y) # geometry grid = CartesianGrid(; origin=(0.0, 0.0), extent=(1.0, 1.0), - size=dims) + dims) free_slip = SBC(0.0, 0.0) free_surf = TBC(0.0, 0.0) @@ -68,8 +69,8 @@ vy(x, y) = -cos(0.5π * x) * sin(0.5π * y) # PT params r = 0.75 re_mech = 2π - lτ_re_m = minimum(extent(grid)) / re_mech - vdτ = minimum(spacing(grid)) / sqrt(ndims(grid)) + lτ_re_m = minimum(extent(grid, Vertex())) / re_mech + vdτ = minimum(spacing(grid, Center(), 1, 1)) / sqrt(ndims(grid)) θ_dτ = lτ_re_m * (r + 4 / 3) / vdτ dτ_r = 1.0 / (θ_dτ + 1.0) nudτ = vdτ * lτ_re_m From 3d7cbdffa30644b64c85dfc8828aa1f09361f753 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 09:43:19 +0100 Subject: [PATCH 03/71] Add doc cleanup --- .github/workflows/DocPreviewCleanup.yml | 34 +++++++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 .github/workflows/DocPreviewCleanup.yml 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 From 93533a8c97bff6a8d699acb76d3f3e962213e969 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 09:48:39 +0100 Subject: [PATCH 04/71] rm exe name --- test/runtests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 32dd4fe2..53a1c09f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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)]...))) From e0812516726b63b4edb07f34d33b6ed673d33cad Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 13:08:46 +0100 Subject: [PATCH 05/71] Fixup tests --- .github/workflows/UnitTests.yml | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/.github/workflows/UnitTests.yml b/.github/workflows/UnitTests.yml index 0cb26030..83b99d99 100644 --- a/.github/workflows/UnitTests.yml +++ b/.github/workflows/UnitTests.yml @@ -59,8 +59,19 @@ 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("https://github.com/PTsolvers/Chmy.jl")' + 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: From 08cb70649822343489f7819b39aaaea271acb2ff Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 13:13:16 +0100 Subject: [PATCH 06/71] Fixup --- .github/workflows/UnitTests.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/UnitTests.yml b/.github/workflows/UnitTests.yml index 83b99d99..1229b9aa 100644 --- a/.github/workflows/UnitTests.yml +++ b/.github/workflows/UnitTests.yml @@ -63,11 +63,11 @@ jobs: # - uses: julia-actions/julia-buildpkg@latest # - uses: julia-actions/julia-runtest@latest # --- testing with un-registered Chmy.jl - - name: run tests + - name: Run tests run: | julia -e 'println("--- :julia: Instantiating project") using Pkg - Pkg.add("https://github.com/PTsolvers/Chmy.jl")' + Pkg.add(url="https://github.com/PTsolvers/Chmy.jl")' julia -e 'println("+++ :julia: Running tests") using Pkg Pkg.test("FastIce"; coverage=true)' From cc7653f3507b9047b471719237a3539c443d7c46 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 13:44:37 +0100 Subject: [PATCH 07/71] Make dummy test pass --- src/FastIce.jl | 2 +- src/LevelSets/compute_level_sets.jl | 2 +- src/LevelSets/signed_distances.jl | 2 +- src/Models/Diffusion/Heat/Heat.jl | 1 + src/Models/Diffusion/Heat/kernels.jl | 4 +- src/Writers.jl | 8 +- test/Project.toml | 1 + test/runtests.jl | 4 +- test/test_chmy_grids.jl | 163 +++++++++++++++++++++++++++ test/test_writers.jl | 5 +- 10 files changed, 176 insertions(+), 16 deletions(-) create mode 100644 test/test_chmy_grids.jl diff --git a/src/FastIce.jl b/src/FastIce.jl index 143c7725..b8408213 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -32,6 +32,6 @@ include("Writers.jl") include("LevelSets/LevelSets.jl") # ice flow models -include("Models/Models.jl") +# include("Models/Models.jl") end # module diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 3bade50b..552dc932 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -11,7 +11,7 @@ end Compute a level set from a given dem. """ -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::CartesianGrid, Ψ_grid::CartesianGrid, +function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::StructuredGrid, Ψ_grid::StructuredGrid, R=LinearAlgebra.I) cutoff = 4maximum(spacing(Ψ_grid)) kernel = _init_level_set!(backend(arch), 256, size(Ψ)) diff --git a/src/LevelSets/signed_distances.jl b/src/LevelSets/signed_distances.jl index 7d4c43f9..5af7c964 100644 --- a/src/LevelSets/signed_distances.jl +++ b/src/LevelSets/signed_distances.jl @@ -51,7 +51,7 @@ 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)) +Base.clamp(p::NTuple{N}, grid::StructuredGrid{N}) where {N} = clamp.(p, Grids.origin(grid), Grids.origin(grid) .+ extent(grid)) function sd_dem(P, cutoff, dem, grid) @inbounds Pp = clamp((P[1], P[2]), grid) diff --git a/src/Models/Diffusion/Heat/Heat.jl b/src/Models/Diffusion/Heat/Heat.jl index 8f97bb8c..b5a0f0dd 100644 --- a/src/Models/Diffusion/Heat/Heat.jl +++ b/src/Models/Diffusion/Heat/Heat.jl @@ -8,6 +8,7 @@ using Chmy.Grids using Chmy.Fields using Chmy.BoundaryConditions using Chmy.KernelLaunch +using Chmy.GridOperators include("kernels.jl") include("boundary_conditions.jl") 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/Writers.jl b/src/Writers.jl index f97204ac..87d6f8c5 100644 --- a/src/Writers.jl +++ b/src/Writers.jl @@ -16,7 +16,7 @@ using MPI Write output `fields` in HDF5 format to a file on `path` for global `grid`. """ -function write_h5(::Architecture, grid::CartesianGrid, path, fields) +function write_h5(::Architecture, grid::StructuredGrid, path, fields) I = CartesianIndices(size(grid)) h5open(path, "w") do io write_dset(io, fields, size(grid), I.indices) @@ -24,7 +24,7 @@ function write_h5(::Architecture, grid::CartesianGrid, path, fields) return end -function write_h5(arch::DistributedArchitecture, grid::CartesianGrid, path, fields) +function write_h5(arch::DistributedArchitecture, grid::StructuredGrid, path, fields) HDF5.has_parallel() || @warn("HDF5 has no parallel support.") topo = topology(arch) comm = cartesian_communicator(topo) @@ -53,7 +53,7 @@ end 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, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) +function write_xdmf(arch::Architecture, grid::StructuredGrid, path, fields, h5_names, timesteps=Float64(0.0)) grid_size = size(grid) grid_spacing = spacing(grid) grid_origin = origin(grid) @@ -64,7 +64,7 @@ function write_xdmf(arch::Architecture, grid::CartesianGrid, path, fields, h5_na return end -function write_xdmf(arch::DistributedArchitecture, grid::CartesianGrid, path, fields, h5_names, timesteps=Float64(0.0)) +function write_xdmf(arch::DistributedArchitecture, grid::StructuredGrid, path, fields, h5_names, timesteps=Float64(0.0)) topo = topology(arch) grid_size = size(grid) grid_spacing = spacing(grid) diff --git a/test/Project.toml b/test/Project.toml index a2881bc9..2650b706 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" diff --git a/test/runtests.jl b/test/runtests.jl index 53a1c09f..08f267d3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,10 +3,10 @@ using FastIce using Pkg -excludedfiles = ["test_excluded.jl"] +excludedfiles = ["test_writers.jl", "test_distributed_writers_3D.jl", "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 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_writers.jl b/test/test_writers.jl index f64d9b50..6727ffac 100644 --- a/test/test_writers.jl +++ b/test/test_writers.jl @@ -32,13 +32,10 @@ XML_ref = """ """ -grid = CartesianGrid(; origin=(-0.4, -0.5, 0.0), - extent=(1.0, 1.1, 1.2), - size=(4, 5, 6)) - for backend in backends @testset "$(basename(@__FILE__)) (backend: $backend)" begin arch = Arch(backend) + grid = UniformGrid(arch; origin=(-0.4, -0.5, 0.0), extent=(1.0, 1.1, 1.2), dims=(4, 5, 6)) Fa = Field(backend, grid, Center()) Fb = Field(backend, grid, Center()) From cb87bd5d5ca58da18f986dfd58ff13c418fa42ea Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 13:47:12 +0100 Subject: [PATCH 08/71] Fixup --- .github/workflows/UnitTests.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/UnitTests.yml b/.github/workflows/UnitTests.yml index 1229b9aa..f61f61cb 100644 --- a/.github/workflows/UnitTests.yml +++ b/.github/workflows/UnitTests.yml @@ -67,7 +67,8 @@ jobs: run: | julia -e 'println("--- :julia: Instantiating project") using Pkg - Pkg.add(url="https://github.com/PTsolvers/Chmy.jl")' + 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)' From 9e25a9f501b64ca6eacc476c88809e36b679641e Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 13:53:21 +0100 Subject: [PATCH 09/71] Add Chmy as not registered --- .buildkite/run_tests.yml | 2 ++ .github/workflows/Documentation.yml | 6 +++++- .github/workflows/UnitTests.yml | 1 + 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/.buildkite/run_tests.yml b/.buildkite/run_tests.yml index 5339ffb4..92b6f00d 100644 --- a/.buildkite/run_tests.yml +++ b/.buildkite/run_tests.yml @@ -13,6 +13,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 +39,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") diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 31fff6f0..43ee5ada 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -18,7 +18,11 @@ jobs: with: version: '1.10' - 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 f61f61cb..c030eb15 100644 --- a/.github/workflows/UnitTests.yml +++ b/.github/workflows/UnitTests.yml @@ -69,6 +69,7 @@ jobs: 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)' From beabc1ebe23e39e58733aa6e42b63e8cbda0be9b Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 4 Mar 2024 13:58:37 +0100 Subject: [PATCH 10/71] Fixup docs --- .buildkite/run_tests.yml | 1 + .github/workflows/Documentation.yml | 1 + .github/workflows/UnitTests.yml | 1 + docs/src/lib/modules.md | 23 +---------------------- 4 files changed, 4 insertions(+), 22 deletions(-) diff --git a/.buildkite/run_tests.yml b/.buildkite/run_tests.yml index 92b6f00d..c94b483e 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: diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 43ee5ada..bf39c935 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: diff --git a/.github/workflows/UnitTests.yml b/.github/workflows/UnitTests.yml index c030eb15..61be3d4d 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: 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 +``` From 626ab8ceb29a8cadfdd211344df3937c816e1e1a Mon Sep 17 00:00:00 2001 From: fquarenghi Date: Tue, 5 Mar 2024 13:59:54 +0100 Subject: [PATCH 11/71] level sets and volume fractions with Chmy backend --- scripts_future_API/Project.toml | 4 +- .../test_levelsets_volumefractions.jl | 109 ++++++++++++++ src/LevelSets/LevelSets.jl | 2 + src/LevelSets/compute_level_sets.jl | 27 ++-- src/LevelSets/compute_volume_fractions.jl | 142 ++++++++++++++++++ src/LevelSets/signed_distances.jl | 21 ++- src/LevelSets/volume_fractions_kernels.jl | 70 +++++++++ 7 files changed, 357 insertions(+), 18 deletions(-) create mode 100644 scripts_future_API/test_levelsets_volumefractions.jl create mode 100644 src/LevelSets/compute_volume_fractions.jl create mode 100644 src/LevelSets/volume_fractions_kernels.jl diff --git a/scripts_future_API/Project.toml b/scripts_future_API/Project.toml index 9b70abea..c6735f55 100644 --- a/scripts_future_API/Project.toml +++ b/scripts_future_API/Project.toml @@ -1,8 +1,8 @@ [deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" 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" diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl new file mode 100644 index 00000000..1de7b789 --- /dev/null +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -0,0 +1,109 @@ +using CUDA +using JLD2 +using FastIce.LevelSets +using Chmy.Grids +using Chmy.Fields +using Chmy.Architectures + + +# Data +# vavilov_path = "../data/vavilov.jld2" +# vavilov_path = "../data/Vavilov/Vavilov_80m.jld2" +vavilov_path = "../data/Vavilov/Vavilov_500m.jld2" +# synthetic_data = "../data/Synthetic/dome_glacier.jld2" +# synthetic_data = "../data/Synthetic/low_res_dome_glacier.jld2" + +# Select backend +backend = CPU() +# backend = CUDABackend() +arch = Arch(backend) + +""" + load_dem_on_GPU(backend, arch::Architecture, path::String) + +Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +initiate grids, copy data on gpu. +""" +function load_dem_on_GPU(backend, arch::Architecture, path::String) + data = load(path) + z_surf = data["DataElevation"].z_surf + z_bed = data["DataElevation"].z_bed + x = data["DataElevation"].x + y = data["DataElevation"].y + nx = length(x) - 1 + ny = length(y) - 1 + nz = 100 + # TODO: choose nz accordingly + surf_lz = maximum(z_surf) - minimum(z_surf) + bed_lz = maximum(z_bed) - minimum(z_bed) + surf_Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(surf_lz, surf_lz, surf_lz), dims=(nx, ny, nz)) + bed_Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_bed)), extent=(bed_lz, bed_lz, bed_lz), dims=(nx, ny, nz)) + surf_dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(surf_lz, surf_lz), dims=(nx, ny)) + bed_dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(bed_lz, bed_lz), dims=(nx, ny)) + surf_field = Field(backend, surf_dem_grid, Vertex()) + bed_field = Field(backend, bed_dem_grid, Vertex()) + set!(surf_field, z_surf) + set!(bed_field, z_bed) + return surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid +end + + +""" + load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) + +Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +initiate grids, copy data on gpu. +""" +function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) + data = load(synthetic_data) + z_surf = data["SyntheticElevation"].z_surf + z_bed = data["SyntheticElevation"].z_bed + nx = size(z_surf)[1] - 1 + ny = size(z_surf)[2] - 1 + nz = 10 + lz = maximum(z_surf) - minimum(z_surf) + Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(lz, lz, lz), dims=(nx, ny, nz)) + # Ψ = Field(backend, Ψ_grid, Vertex()) + dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(lz, lz), dims=(nx, ny)) + dem_bed = Field(backend, dem_grid, Vertex()) + dem_surf = Field(backend, dem_grid, Vertex()) + set!(dem_bed, z_bed) + set!(dem_surf, z_surf) + return dem_surf, dem_bed, dem_grid, Ψ_grid +end + +# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem_on_GPU(backend, arch, synthetic_data); +surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem_on_GPU(backend, arch, vavilov_path); + +Ψ = ( + na=Field(backend, surf_Ψ_grid, Vertex()), + ns=Field(backend, bed_Ψ_grid, Vertex()), +) + +compute_level_set_from_dem!(backend, Ψ[1], surf_field, surf_dem_grid, surf_Ψ_grid) +compute_level_set_from_dem!(backend, Ψ[2], bed_field, bed_dem_grid, bed_Ψ_grid) + +# for phase in eachindex(Ψ) +# compute_level_set_from_dem!(backend, Ψ[phase], dem_surf, dem_grid, Ψ_grid) +# end + +wt = ( + na=volfrac_field(backend, surf_Ψ_grid), + ns=volfrac_field(backend, bed_Ψ_grid), +) + +compute_volume_fractions_from_level_set!(backend, wt[1], Ψ[1], surf_Ψ_grid) +compute_volume_fractions_from_level_set!(backend, wt[2], Ψ[2], bed_Ψ_grid) + + +# for phase in eachindex(Ψ) +# compute_volume_fractions_from_level_set!(backend, wt[phase], Ψ[phase], Ψ_grid) +# end + +# Save +jldopen(vavilov_path, "a+") do file + file["level_sets_na"] = Array(interior(Ψ[1])) + file["level_sets_ns"] = Array(interior(Ψ[2])) + file["volume_frac_na"] = Array.(interior.(Tuple(wt[1]))) + file["volume_frac_ns"] = Array.(interior.(Tuple(wt[2]))) +end diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index 5eb8428c..37ad3edf 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -1,6 +1,7 @@ module LevelSets export compute_level_set_from_dem! +export volfrac_field, compute_volume_fractions_from_level_set! using Chmy using Chmy.Grids @@ -11,5 +12,6 @@ using LinearAlgebra, GeometryBasics include("signed_distances.jl") include("compute_level_sets.jl") +include("compute_volume_fractions.jl") end diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 552dc932..3b6b6b36 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -1,20 +1,25 @@ -@kernel inbounds = true function _init_level_set!(Ψ::Field, dem::Field, dem_grid::StructuredGrid, Ψ_grid::StructuredGrid, cutoff, R) - I = @index(Global, Cartesian) - x, y, z = coord(Ψ_grid, location(Ψ), I) +""" + _init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff::AbstractFloat, R::AbstractMatrix) + +Initialize level sets. +""" +@kernel function init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) + I = @index(Global, NTuple) + 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 + Ψ[I...] = ud * sgn end + """ - compute_level_set_from_dem!(arch, Ψ, dem, dem_grid, Ψ_grid, R[=LinearAlgebra.I]) + compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) -Compute a level set from a given dem. +Compute level sets from dem. """ -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::StructuredGrid, Ψ_grid::StructuredGrid, - R=LinearAlgebra.I) - cutoff = 4maximum(spacing(Ψ_grid)) - kernel = _init_level_set!(backend(arch), 256, size(Ψ)) +function compute_level_set_from_dem!(backend, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) + kernel = init_level_set!(backend, 256, size(Ψ)) + cutoff = 4maximum(spacing(Ψ_grid, Center(), 1, 1, 1)) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) return -end +end \ No newline at end of file diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl new file mode 100644 index 00000000..37c8dbb5 --- /dev/null +++ b/src/LevelSets/compute_volume_fractions.jl @@ -0,0 +1,142 @@ +export volfrac + +@inline perturb(ϕ) = abs(ϕ) > 1e-20 ? ϕ : (ϕ > 0 ? 1e-20 : -1e-20) + +@inline trivol(v1, v2, v3) = 0.5 * abs(cross(v3 - v1, v2 - v1)) + +function volfrac(tri, ϕ::Vec3{T})::T where {T} + v1, v2, v3 = tri + if ϕ[1] < 0 && ϕ[2] < 0 && ϕ[3] < 0 # --- + return trivol(v1, v2, v3) + elseif ϕ[1] > 0 && ϕ[2] > 0 && ϕ[3] > 0 # +++ + return 0.0 + end + @inline vij(i, j) = tri[j] * (ϕ[i] / (ϕ[i] - ϕ[j])) - tri[i] * (ϕ[j] / (ϕ[i] - ϕ[j])) + v12, v13, v23 = vij(1, 2), vij(1, 3), vij(2, 3) + if ϕ[1] < 0 + if ϕ[2] < 0 + trivol(v1, v23, v13) + trivol(v1, v2, v23) # --+ + else + if ϕ[3] < 0 + trivol(v3, v12, v23) + trivol(v3, v1, v12) # -+- + else + trivol(v1, v12, v13) # -++ + end + end + else + if ϕ[2] < 0 + if ϕ[3] < 0 + trivol(v2, v13, v12) + trivol(v2, v3, v13) # +-- + else + trivol(v12, v2, v23) # +-+ + end + else + trivol(v13, v23, v3) # ++- + end + end +end + +function volfrac(rect::Rect2{T}, ϕ::Vec4{T}) where {T} + or, ws = Grids.origin(rect), widths(rect) + v1, v2, v3, v4 = or, or + Vec(ws[1], 0.0), or + ws, or + Vec(0.0, ws[2]) + ϕ1, ϕ2, ϕ3, ϕ4 = perturb.(ϕ) + return volfrac(Vec(v1, v2, v3), Vec3{T}(ϕ1, ϕ2, ϕ3)) + + volfrac(Vec(v1, v3, v4), Vec3{T}(ϕ1, ϕ3, ϕ4)) +end + +@inline tetvol(v1, v2, v3, v4) = abs(det([v2 - v1 v3 - v1 v4 - v1])) / 6.0 + +function volfrac(tet, ϕ::Vec4) + v1, v2, v3, v4 = tet + @inline vij(i, j) = tet[j] * (ϕ[i] / (ϕ[i] - ϕ[j])) - tet[i] * (ϕ[j] / (ϕ[i] - ϕ[j])) + nneg = count(ϕ .< 0) + if nneg == 0 # ++++ + return 0.0 + elseif nneg == 1 # -+++ + if ϕ[1] < 0 + return tetvol(v1, vij(1, 2), vij(1, 3), vij(1, 4)) + elseif ϕ[2] < 0 + return tetvol(v2, vij(2, 1), vij(2, 3), vij(2, 4)) + elseif ϕ[3] < 0 + return tetvol(v3, vij(3, 1), vij(3, 2), vij(3, 4)) + else # ϕ[4] < 0 + return tetvol(v4, vij(4, 1), vij(4, 2), vij(4, 3)) + end + elseif nneg == 2 # --++ + if ϕ[1] < 0 && ϕ[2] < 0 + return tetvol(v1, v2, vij(1, 3), vij(2, 4)) + + tetvol(vij(2, 3), v2, vij(1, 3), vij(2, 4)) + + tetvol(v1, vij(1, 4), vij(1, 3), vij(2, 4)) + elseif ϕ[1] < 0 && ϕ[3] < 0 + return tetvol(v1, v3, vij(1, 4), vij(3, 2)) + + tetvol(vij(3, 4), v3, vij(1, 4), vij(3, 2)) + + tetvol(v1, vij(1, 2), vij(1, 4), vij(3, 2)) + elseif ϕ[1] < 0 && ϕ[4] < 0 + return tetvol(v1, v4, vij(1, 2), vij(4, 3)) + + tetvol(vij(4, 2), v4, vij(1, 2), vij(4, 3)) + + tetvol(v1, vij(1, 3), vij(1, 2), vij(4, 3)) + elseif ϕ[2] < 0 && ϕ[3] < 0 + return tetvol(v3, v2, vij(3, 1), vij(2, 4)) + + tetvol(vij(2, 1), v2, vij(3, 1), vij(2, 4)) + + tetvol(v3, vij(3, 4), vij(3, 1), vij(2, 4)) + elseif ϕ[2] < 0 && ϕ[4] < 0 + return tetvol(v4, v2, vij(4, 1), vij(2, 3)) + + tetvol(vij(2, 1), v2, vij(4, 1), vij(2, 3)) + + tetvol(v4, vij(4, 3), vij(4, 1), vij(2, 3)) + else # ϕ[3] < 0 && ϕ[4] < 0 + return tetvol(v3, v4, vij(3, 1), vij(4, 2)) + + tetvol(vij(4, 1), v4, vij(3, 1), vij(4, 2)) + + tetvol(v3, vij(3, 2), vij(3, 1), vij(4, 2)) + end + elseif nneg == 3 # ---+ + vol_tot = tetvol(v1, v2, v3, v4) + if ϕ[1] >= 0 + return vol_tot - tetvol(v1, vij(1, 2), vij(1, 3), vij(1, 4)) + elseif ϕ[2] >= 0 + return vol_tot - tetvol(v2, vij(2, 1), vij(2, 3), vij(2, 4)) + elseif ϕ[3] >= 0 + return vol_tot - tetvol(v3, vij(3, 1), vij(3, 2), vij(3, 4)) + else # ϕ[4] >= 0 + return vol_tot - tetvol(v4, vij(4, 1), vij(4, 2), vij(4, 3)) + end + else # ---- + return tetvol(v1, v2, v3, v4) + end +end + +function volfrac(rect::Rect3, ϕ::Vec{8}) + or, ws = GeometryBasics.origin(rect), widths(rect) + v000, v001, v100, v101 = or, or + Vec(ws[1], 0.0, 0.0), or + Vec(0.0, ws[2], 0.0), or + Vec(ws[1], ws[2], 0.0) + v010, v011, v110, v111 = or + Vec(0.0, 0.0, ws[3]), or + Vec(ws[1], 0.0, ws[3]), or + Vec(0.0, ws[2], ws[3]), + or + Vec(ws[1], ws[2], ws[3]) + ϕ = perturb.(ϕ) + return volfrac(Vec(v000, v100, v010, v001), Vec(ϕ[1], ϕ[5], ϕ[3], ϕ[2])) + + volfrac(Vec(v110, v100, v010, v111), Vec(ϕ[7], ϕ[5], ϕ[3], ϕ[7])) + + volfrac(Vec(v101, v100, v111, v001), Vec(ϕ[6], ϕ[5], ϕ[7], ϕ[2])) + + volfrac(Vec(v011, v111, v010, v001), Vec(ϕ[4], ϕ[7], ϕ[3], ϕ[2])) + + volfrac(Vec(v111, v100, v010, v001), Vec(ϕ[7], ϕ[5], ϕ[3], ϕ[2])) +end + +function volfrac_field(backend::Backend, grid::UniformGrid{2}, args...; kwargs...) + (c=Field(backend, grid, Center()), + x=Field(backend, grid, (Vertex(), Center())), + y=Field(backend, grid, (Center(), Vertex())), + xy=Field(backend, grid, Vertex())) +end + +function volfrac_field(backend::Backend, grid::UniformGrid{3}, args...; kwargs...) + (c=Field(backend, grid, Center()), + x=Field(backend, grid, (Vertex(), Center(), Center())), + y=Field(backend, grid, (Center(), Vertex(), Center())), + z=Field(backend, grid, (Center(), Center(), Vertex())), + xy=Field(backend, grid, (Vertex(), Vertex(), Center())), + xz=Field(backend, grid, (Vertex(), Center(), Vertex())), + yz=Field(backend, grid, (Center(), Vertex(), Vertex()))) +end + +include("volume_fractions_kernels.jl") + +function compute_volume_fractions_from_level_set!(backend::Backend, wt, Ψ, grid::UniformGrid) + kernel = kernel_compute_volume_fractions_from_level_set!(backend, 256, size(Ψ)) + kernel(wt, Ψ, grid) +end \ No newline at end of file diff --git a/src/LevelSets/signed_distances.jl b/src/LevelSets/signed_distances.jl index 5af7c964..36d97ee3 100644 --- a/src/LevelSets/signed_distances.jl +++ b/src/LevelSets/signed_distances.jl @@ -1,5 +1,8 @@ +export sd_dem + @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) @@ -23,18 +26,22 @@ 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 end + @inline inc(I, dim) = Base.setindex(I, I[dim] + 1, dim) @inline inc(I) = I + oneunit(I) + @inline function triangle_pair(Iv, dem, grid) @inline function sample_dem(I) @inbounds x, y = coord(grid, location(dem), I) @@ -45,13 +52,16 @@ end return T_BL, T_TR end + @inline function distance_to_triangle_pair(P, Iv, dem, grid) T_BL, T_TR = triangle_pair(Iv, dem, grid) ud = min(ud_triangle(P, T_BL...), ud_triangle(P, T_TR...)) return ud, sign_triangle(P, T_BL...) end -Base.clamp(p::NTuple{N}, grid::StructuredGrid{N}) where {N} = clamp.(p, Grids.origin(grid), Grids.origin(grid) .+ extent(grid)) + +Base.clamp(p::NTuple{N}, grid::UniformGrid{N}) where {N} = clamp.(p, Grids.origin(grid, Vertex()), Grids.origin(grid, Vertex()) .+ extent(grid, Vertex())) + function sd_dem(P, cutoff, dem, grid) @inbounds Pp = clamp((P[1], P[2]), grid) @@ -68,3 +78,4 @@ function sd_dem(P, cutoff, dem, grid) end return ud, sgn end + diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl new file mode 100644 index 00000000..dd16cbf9 --- /dev/null +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -0,0 +1,70 @@ +# Volume fraction kernel (2D) +@kernel function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid) + ix, iy = @index(Global, NTuple) + dx, dy = spacing(grid, Center(), 1, 1) + cell = Rect(Vec(0.0, 0.0), Vec(dx, dy)) + ω = GeometryBasics.volume(cell) + Ψ_ax(dix, diy) = 0.5 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix+1, iy+diy]) + Ψ_ay(dix, diy) = 0.5 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix, iy+diy+1]) + Ψ_axy(dix, diy) = 0.25 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix+1, iy+diy] + Ψ[ix+dix, iy+diy+1] + Ψ[ix+dix+1, iy+diy+1]) + # cell centers + Ψs = Vec{4}(Ψ[ix, iy], Ψ[ix+1, iy], Ψ[ix+1, iy+1], Ψ[ix, iy+1]) + wt.c[ix, iy] = volfrac(cell, Ψs) / ω + # x faces + Ψs = Vec{4}(Ψ_ax(0, 0), Ψ_ax(1, 0), Ψ_ax(1, 1), Ψ_ax(0, 1)) + wt.x[ix, iy] = volfrac(cell, Ψs) / ω + # y faces + Ψs = Vec{4}(Ψ_ay(0, 0), Ψ_ay(1, 0), Ψ_ay(1, 1), Ψ_ay(0, 1)) + wt.y[ix, iy] = volfrac(cell, Ψs) / ω + # xy edges + Ψs = Vec{4}(Ψ_axy(0, 0), Ψ_axy(1, 0), Ψ_axy(1, 1), Ψ_axy(0, 1)) + wt.xy[ix, iy] = volfrac(cell, Ψs) / ω +end + +# Volume fraction kernel (3D) +@kernel function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid) + ix, iy, iz = @index(Global, NTuple) + dx, dy, dz = spacing(grid, Center(), 1, 1, 1) + cell = Rect(Vec(0.0, 0.0, 0.0), Vec(dx, dy, dz)) + ω = GeometryBasics.volume(cell) + Ψ_ax(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz]) + Ψ_ay(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz]) + Ψ_az(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix, iy+diy, iz+diz+1]) + function Ψ_axy(dix, diy, diz) + 0.25 * + (Ψ[ix+dix, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + end + function Ψ_axz(dix, diy, diz) + 0.25 * + (Ψ[ix+dix, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + end + function Ψ_ayz(dix, diy, diz) + 0.25 * + (Ψ[ix+dix+1, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + end + # cell centers + Ψs = Vec{8}(Ψ[ix, iy, iz], Ψ[ix+1, iy, iz], Ψ[ix, iy+1, iz], Ψ[ix+1, iy+1, iz], Ψ[ix, iy, iz+1], Ψ[ix+1, iy, iz+1], Ψ[ix, iy+1, iz+1], + Ψ[ix+1, iy+1, iz+1]) + wt.c[ix, iy, iz] = volfrac(cell, Ψs) / ω + # x faces + Ψs = Vec{8}(Ψ_ax(0, 0, 0), Ψ_ax(1, 0, 0), Ψ_ax(0, 1, 0), Ψ_ax(1, 1, 0), Ψ_ax(0, 0, 1), Ψ_ax(1, 0, 1), Ψ_ax(0, 1, 1), Ψ_ax(1, 1, 1)) + wt.x[ix, iy, iz] = volfrac(cell, Ψs) / ω + # y faces + Ψs = Vec{8}(Ψ_ay(0, 0, 0), Ψ_ay(1, 0, 0), Ψ_ay(0, 1, 0), Ψ_ay(1, 1, 0), Ψ_ay(0, 0, 1), Ψ_ay(1, 0, 1), Ψ_ay(0, 1, 1), Ψ_ay(1, 1, 1)) + wt.y[ix, iy, iz] = volfrac(cell, Ψs) / ω + # z faces + Ψs = Vec{8}(Ψ_az(0, 0, 0), Ψ_az(1, 0, 0), Ψ_az(0, 1, 0), Ψ_az(1, 1, 0), Ψ_az(0, 0, 1), Ψ_az(1, 0, 1), Ψ_az(0, 1, 1), Ψ_az(1, 1, 1)) + wt.z[ix, iy, iz] = volfrac(cell, Ψs) / ω + # xy edges + Ψs = Vec{8}(Ψ_axy(0, 0, 0), Ψ_axy(1, 0, 0), Ψ_axy(0, 1, 0), Ψ_axy(1, 1, 0), Ψ_axy(0, 0, 1), Ψ_axy(1, 0, 1), Ψ_axy(0, 1, 1), + Ψ_axy(1, 1, 1)) + wt.xy[ix, iy, iz] = volfrac(cell, Ψs) / ω + # xz edges + Ψs = Vec{8}(Ψ_axz(0, 0, 0), Ψ_axz(1, 0, 0), Ψ_axz(0, 1, 0), Ψ_axz(1, 1, 0), Ψ_axz(0, 0, 1), Ψ_axz(1, 0, 1), Ψ_axz(0, 1, 1), + Ψ_axz(1, 1, 1)) + wt.xz[ix, iy, iz] = volfrac(cell, Ψs) / ω + # yz edges + Ψs = Vec{8}(Ψ_ayz(0, 0, 0), Ψ_ayz(1, 0, 0), Ψ_ayz(0, 1, 0), Ψ_ayz(1, 1, 0), Ψ_ayz(0, 0, 1), Ψ_ayz(1, 0, 1), Ψ_ayz(0, 1, 1), + Ψ_ayz(1, 1, 1)) + wt.yz[ix, iy, iz] = volfrac(cell, Ψs) / ω +end From ca353e008969df51d67e0cec4255b2c219a0d8ef Mon Sep 17 00:00:00 2001 From: Filippo Quarenghi Date: Tue, 5 Mar 2024 16:28:06 +0100 Subject: [PATCH 12/71] Update src/LevelSets/compute_level_sets.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Ludovic Räss <61313342+luraess@users.noreply.github.com> --- src/LevelSets/compute_level_sets.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 3b6b6b36..e486cedd 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -19,7 +19,7 @@ Compute level sets from dem. """ function compute_level_set_from_dem!(backend, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) kernel = init_level_set!(backend, 256, size(Ψ)) - cutoff = 4maximum(spacing(Ψ_grid, Center(), 1, 1, 1)) + cutoff = 4maximum(spacing(Ψ_grid, Center())) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) return end \ No newline at end of file From 766136bf215156698c37700de2a81268e708685a Mon Sep 17 00:00:00 2001 From: fquarenghi Date: Tue, 5 Mar 2024 16:31:45 +0100 Subject: [PATCH 13/71] fixed names and toml file --- scripts_future_API/Project.toml | 3 --- .../test_levelsets_volumefractions.jl | 16 +++++++--------- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/scripts_future_API/Project.toml b/scripts_future_API/Project.toml index c6735f55..94563595 100644 --- a/scripts_future_API/Project.toml +++ b/scripts_future_API/Project.toml @@ -1,9 +1,6 @@ [deps] -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" FastIce = "e0de9f13-a007-490e-b696-b07d031015ca" -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/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index 1de7b789..965fa1a7 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -4,7 +4,7 @@ using FastIce.LevelSets using Chmy.Grids using Chmy.Fields using Chmy.Architectures - +using KernelAbstractions # Data # vavilov_path = "../data/vavilov.jld2" @@ -19,12 +19,12 @@ backend = CPU() arch = Arch(backend) """ - load_dem_on_GPU(backend, arch::Architecture, path::String) + load_dem(backend, arch::Architecture, path::String) Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, initiate grids, copy data on gpu. """ -function load_dem_on_GPU(backend, arch::Architecture, path::String) +function load_dem(backend, arch::Architecture, path::String) data = load(path) z_surf = data["DataElevation"].z_surf z_bed = data["DataElevation"].z_bed @@ -47,14 +47,13 @@ function load_dem_on_GPU(backend, arch::Architecture, path::String) return surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid end - """ - load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) + load_synth_dem(backend, arch::Architecture, synthetic_data::String) Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, initiate grids, copy data on gpu. """ -function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) +function load_synth_dem(backend, arch::Architecture, synthetic_data::String) data = load(synthetic_data) z_surf = data["SyntheticElevation"].z_surf z_bed = data["SyntheticElevation"].z_bed @@ -63,7 +62,6 @@ function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::Stri nz = 10 lz = maximum(z_surf) - minimum(z_surf) Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(lz, lz, lz), dims=(nx, ny, nz)) - # Ψ = Field(backend, Ψ_grid, Vertex()) dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(lz, lz), dims=(nx, ny)) dem_bed = Field(backend, dem_grid, Vertex()) dem_surf = Field(backend, dem_grid, Vertex()) @@ -72,8 +70,8 @@ function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::Stri return dem_surf, dem_bed, dem_grid, Ψ_grid end -# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem_on_GPU(backend, arch, synthetic_data); -surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem_on_GPU(backend, arch, vavilov_path); +# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem(backend, arch, synthetic_data); +surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem(backend, arch, vavilov_path); Ψ = ( na=Field(backend, surf_Ψ_grid, Vertex()), From c94b240c32c5761316903e188762c00bcaf7204c Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 5 Mar 2024 17:06:30 +0100 Subject: [PATCH 14/71] Fix writers --- src/Writers.jl | 35 ++++++------ test/runtests.jl | 2 +- test/test_distributed_writers_3D.jl | 88 +++++++++++++---------------- test/test_writers.jl | 45 +++++++-------- 4 files changed, 83 insertions(+), 87 deletions(-) diff --git a/src/Writers.jl b/src/Writers.jl index 87d6f8c5..900823e2 100644 --- a/src/Writers.jl +++ b/src/Writers.jl @@ -16,25 +16,26 @@ using MPI Write output `fields` in HDF5 format to a file on `path` for global `grid`. """ -function write_h5(::Architecture, grid::StructuredGrid, path, fields) - I = CartesianIndices(size(grid)) +function write_h5(::Architecture, grid::UniformGrid, path, fields) + I = CartesianIndices(size(grid, Center())) h5open(path, "w") do io - write_dset(io, fields, size(grid), I.indices) + write_dset(io, fields, size(grid, Center()), I.indices) end return end -function write_h5(arch::DistributedArchitecture, grid::StructuredGrid, path, fields) +function write_h5(arch::DistributedArchitecture, grid::UniformGrid, path, fields) HDF5.has_parallel() || @warn("HDF5 has no parallel support.") topo = topology(arch) - comm = cartesian_communicator(topo) - coords = coordinates(topo) - sz = size(local_grid(grid, topo)) + comm = cart_comm(topo) + coords = Distributed.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) + write_dset(io, fields, global_grid_size, I.indices) end return end @@ -53,10 +54,10 @@ end 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, grid::StructuredGrid, path, fields, h5_names, timesteps=Float64(0.0)) - 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(grid) + grid_origin = origin(grid, Center()) xdoc = generate_xdmf(grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) @@ -64,13 +65,13 @@ function write_xdmf(arch::Architecture, grid::StructuredGrid, path, fields, h5_n return end -function write_xdmf(arch::DistributedArchitecture, grid::StructuredGrid, path, fields, h5_names, timesteps=Float64(0.0)) +function write_xdmf(arch::DistributedArchitecture, grid::UniformGrid, path, fields, h5_names, timesteps=Float64(0.0)) topo = topology(arch) - grid_size = size(grid) + global_grid_size = dims(topo) .* 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) + xdoc = generate_xdmf(global_grid_size, grid_spacing, grid_origin, fields, h5_names, timesteps) save_file(xdoc, path) return @@ -102,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") @@ -112,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/runtests.jl b/test/runtests.jl index 08f267d3..11e776c5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using FastIce using Pkg -excludedfiles = ["test_writers.jl", "test_distributed_writers_3D.jl", "test_excluded.jl"] # tmp as WIP +excludedfiles = ["test_excluded.jl"] # tmp as WIP # distributed test_distributed = ["test_distributed_writers_3D.jl"] 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_writers.jl b/test/test_writers.jl index 6727ffac..904ee3a0 100644 --- a/test/test_writers.jl +++ b/test/test_writers.jl @@ -3,6 +3,7 @@ include("common.jl") using Chmy.Architectures using Chmy.Fields using Chmy.Grids + using FastIce.Writers using HDF5 @@ -17,14 +18,14 @@ XML_ref = """