Update routines for la mem

[IskanderI]

No description provided.

[github-actions]

Your PR requires formatting changes to meet the project's style guidelines.
Please consider running Runic (git runic main) to apply these changes.

Click here to view the suggested changes.

diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl
index 3095f7e..3e8cecc 100644
--- a/src/LaMEM_io.jl
+++ b/src/LaMEM_io.jl
@@ -37,7 +37,7 @@ struct LaMEMPartitioningInfo <: AbstractGeneralGrid
     ix::Vector{Int64}
     iy::Vector{Int64}
     iz::Vector{Int64}
-    
+
 
 end
 
@@ -1054,14 +1054,14 @@ The file is named `ProcessorPartitioning_(P.nProcX*P.nProcY*P.nProcZ)cpu_(P.nPro
 The coordinates are written as Float64, and the processor counts and node counts as Int64 or Int32, depending on the `is64bit` flag.
 """
 function write_processor_partitioning_LaMEM(
-    P::LaMEMPartitioningInfo;
-    is64bit::Bool = false
-)
-    xcoor = range(P.xc[1], P.xc[end], length=P.nNodeX)
-    ycoor = range(P.yc[1], P.yc[end], length=P.nNodeY)
-    zcoor = range(P.zc[1], P.zc[end], length=P.nNodeZ)
-
-    filename = "ProcessorPartitioning_$(P.nProcX*P.nProcY*P.nProcZ)cpu_$(P.nProcX).$(P.nProcY).$(P.nProcZ).bin"
+        P::LaMEMPartitioningInfo;
+        is64bit::Bool = false
+    )
+    xcoor = range(P.xc[1], P.xc[end], length = P.nNodeX)
+    ycoor = range(P.yc[1], P.yc[end], length = P.nNodeY)
+    zcoor = range(P.zc[1], P.zc[end], length = P.nNodeZ)
+
+    filename = "ProcessorPartitioning_$(P.nProcX * P.nProcY * P.nProcZ)cpu_$(P.nProcX).$(P.nProcY).$(P.nProcZ).bin"
     typ = is64bit ? Int64 : Int32
     open(filename, "w") do io
         # Write processor counts (Int64, big-endian)
@@ -1164,13 +1164,15 @@ function Base.show(io::IO, d::LaMEMPartitioningInfo)
     println(io, "  ix     : $(d.ix)")
     println(io, "  iy     : $(d.iy)")
     println(io, "  iz     : $(d.iz)")
-    
+
     return nothing
 
 end
 
 function check_markers_directory(directory)
-    return if !isdir(directory) mkdir(directory) end
+    return if !isdir(directory)
+        mkdir(directory)
+    end
 end
 
 function get_LaMEM_grid_info(file; args::Union{String, Nothing} = nothing)
@@ -1304,18 +1306,18 @@ function get_proc_grid(Grid_info, p_dist, proc_bounds, proc_num, RandomNoise)
         dXNoise = zeros(size(X)) + dx
         dYNoise = zeros(size(Y)) + dy
         dZNoise = zeros(size(Z)) + dz
-        
+
         dXNoise = dXNoise .* (rand(size(dXNoise)) - 0.5)
         dYNoise = dYNoise .* (rand(size(dYNoise)) - 0.5)
         dZNoise = dZNoise .* (rand(size(dZNoise)) - 0.5)
-        
+
         X .+= dXNoise
         Y .+= dYNoise
         Z .+= dZNoise
         x = X(1, :, 1)
         y = Y(:, 1, 1)
         z = Z(1, 1, :)
-    
+
     end
 
     Grid = LaMEM_grid(
@@ -1467,54 +1469,54 @@ function decompose_mpi_ranks(total_ranks::Int, nx::Int, ny::Int, nz::Int)
     for px in factors
         remaining = total_ranks ÷ px
         rem_factors = get_factors(remaining)
-        
+
         for py in rem_factors
             pz = remaining ÷ py
-            
+
             # Skip invalid combinations
             if px * py * pz != total_ranks
                 continue
             end
-            
+
             # Calculate local grid sizes
             local_nx = nx / px
             local_ny = ny / py
             local_nz = nz / pz
-            
+
             if mod(nx, px) != 0 || mod(ny, py) != 0 || mod(nz, pz) != 0
                 continue  # Skip this iteration if any of them is not an integer
             end
-    
+
             # Calculate aspect ratios (always ≥ 1.0)
-            ar_xy = max(local_nx/local_ny, local_ny/local_nx)
-            ar_xz = max(local_nx/local_nz, local_nz/local_nx)
-            ar_yz = max(local_ny/local_nz, local_nz/local_ny)
-            
+            ar_xy = max(local_nx / local_ny, local_ny / local_nx)
+            ar_xz = max(local_nx / local_nz, local_nz / local_nx)
+            ar_yz = max(local_ny / local_nz, local_nz / local_ny)
+
             # Metric: average deviation from aspect ratio of 1.0
             metric = (ar_xy + ar_xz + ar_yz) / 3.0
 
             # Update best configuration if this one is better
-        if metric < best_metric
+            if metric < best_metric
 
-            best_metric = metric
-            best_px = px
-            best_py = py
-            best_pz = pz
-            first_best_metric = metric  # Store the first best metric
+                best_metric = metric
+                best_px = px
+                best_py = py
+                best_pz = pz
+                first_best_metric = metric  # Store the first best metric
 
-        elseif isapprox(metric, best_metric, rtol=1e-10)
+            elseif isapprox(metric, best_metric, rtol = 1.0e-10)
 
-            # If the metric is equal to the first best metric, skip updating
-            if first_best_metric !== nothing && isapprox(metric, first_best_metric, rtol=1e-10)
-                continue
-            end
+                # If the metric is equal to the first best metric, skip updating
+                if first_best_metric !== nothing && isapprox(metric, first_best_metric, rtol = 1.0e-10)
+                    continue
+                end
 
-        end
+            end
         end
     end
     if best_metric == Inf
         error("No valid decomposition found for total_ranks = $total_ranks")
-        
+
     end
     return (best_px, best_py, best_pz)
 end
@@ -1527,8 +1529,8 @@ function get_factors(n::Int)
     for i in 1:isqrt(n)
         if n % i == 0
             push!(factors, i)
-            if i != n÷i
-                push!(factors, n÷i)
+            if i != n ÷ i
+                push!(factors, n ÷ i)
             end
         end
     end
@@ -1553,47 +1555,49 @@ Parameters:
 Returns:
 - ModelDomain struct containing all domain decomposition information
 """
-function setup_model_domain(coord_x::AbstractVector{<:Real}, 
-                            coord_y::AbstractVector{<:Real}, 
-                            coord_z::AbstractVector{<:Real},
-                            nx::Int, ny::Int, nz::Int,
-                            n_ranks::Int; verbose::Bool = true)
-    
+function setup_model_domain(
+        coord_x::AbstractVector{<:Real},
+        coord_y::AbstractVector{<:Real},
+        coord_z::AbstractVector{<:Real},
+        nx::Int, ny::Int, nz::Int,
+        n_ranks::Int; verbose::Bool = true
+    )
+
     # Verify input vectors have correct size
     if any(length.([coord_x, coord_y, coord_z]) .!= 2)
         error("coord_x, coord_y, and coord_z must be 2-element vectors [min, max]")
     end
-    
+
     # Generate full coordinate vectors
     nnodx = nx + 1
     nnody = ny + 1
     nnodz = nz + 1
-    
-    xcoor = range(coord_x[1], coord_x[2], length=nnodx)
-    ycoor = range(coord_y[1], coord_y[2], length=nnody)
-    zcoor = range(coord_z[1], coord_z[2], length=nnodz)
-    
+
+    xcoor = range(coord_x[1], coord_x[2], length = nnodx)
+    ycoor = range(coord_y[1], coord_y[2], length = nnody)
+    zcoor = range(coord_z[1], coord_z[2], length = nnodz)
+
     check_multigrid_compatibility(nx, ny, nz, n_ranks, verbose)
 
     # Decompose MPI ranks into 3D processor grid
     Nprocx, Nprocy, Nprocz = decompose_mpi_ranks(n_ranks, nx, ny, nz)
-    
+
     # Calculate subdomain divisions
     function calculate_domain_divisions(N::Int, nproc::Int)
         base_size = div(N, nproc)
         remainder = N % nproc
-        
+
         indices = zeros(Int, nproc + 1)
         indices[1] = 1
-        
+
         for i in 1:nproc
             local_size = base_size + (i <= remainder ? 1 : 0)
             indices[i + 1] = indices[i] + local_size
         end
-        
+
         return indices
     end
-    
+
     # Calculate divisions for each direction
     ix = calculate_domain_divisions(nx, Nprocx)
     iy = calculate_domain_divisions(ny, Nprocy)
@@ -1601,14 +1605,14 @@ function setup_model_domain(coord_x::AbstractVector{<:Real},
 
     P = LaMEMPartitioningInfo(
         Nprocx, Nprocy, Nprocz,
-        nnodx, nnody, nnodz, 
+        nnodx, nnody, nnodz,
         xcoor[ix], ycoor[iy], zcoor[iz],
         ix, iy, iz
-        )
+    )
 
-        xcoor = []
-        ycoor = []
-        zcoor = []
+    xcoor = []
+    ycoor = []
+    zcoor = []
 
     return P
 
@@ -1623,9 +1627,9 @@ Check if the model resolution and MPI ranks are optimal for multigrid methods.
 function check_multigrid_compatibility(nx::Int, ny::Int, nz::Int, total_ranks::Int, verbose::Bool = true)
     # Generate arrays of valid numbers for resolution
     powers_of_2 = [2^i for i in 1:15]  # Up to 32768
-    multiples_of_48 = [48*i for i in 1:683]  # Up to 32784
+    multiples_of_48 = [48 * i for i in 1:683]  # Up to 32784
     valid_resolution_numbers = sort(unique([powers_of_2; multiples_of_48]))
-    
+
     # Generate valid rank sequence algorithmically
     valid_rank_numbers = generate_valid_rank_sequence()
 
@@ -1637,12 +1641,12 @@ function check_multigrid_compatibility(nx::Int, ny::Int, nz::Int, total_ranks::I
         elseif idx == 1
             return valid_array[1]
         else
-            prev = valid_array[idx-1]
+            prev = valid_array[idx - 1]
             curr = valid_array[idx]
             return abs(n - prev) < abs(n - curr) ? prev : curr
         end
     end
-    
+
     # Function to check if a number is valid for resolution
     function is_valid_resolution(n)
         # Check if power of 2
@@ -1652,20 +1656,20 @@ function check_multigrid_compatibility(nx::Int, ny::Int, nz::Int, total_ranks::I
         # Check if multiple of 48
         return n % 48 == 0
     end
-    
+
     # Function to check if a number is valid for ranks
     function is_valid_ranks(n)
         return n in valid_rank_numbers
     end
-    
+
     # Check each value
     is_nx_valid = is_valid_resolution(nx)
     is_ny_valid = is_valid_resolution(ny)
     is_nz_valid = is_valid_resolution(nz)
     is_ranks_valid = is_valid_ranks(total_ranks)
-    
+
     all_valid = is_nx_valid && is_ny_valid && is_nz_valid && is_ranks_valid
-    
+
     # Generate suggestions for invalid values
     suggestions = Dict{String, Int}()
     if !is_nx_valid
@@ -1680,7 +1684,7 @@ function check_multigrid_compatibility(nx::Int, ny::Int, nz::Int, total_ranks::I
     if !is_ranks_valid
         suggestions["total_ranks"] = find_closest_valid(total_ranks, valid_rank_numbers)
     end
-    
+
     # Print warnings and suggestions
     if !all_valid && verbose
         println("\nWARNING: Non-optimal configuration detected for multigrid method!")
@@ -1689,14 +1693,14 @@ function check_multigrid_compatibility(nx::Int, ny::Int, nz::Int, total_ranks::I
         println("ny = $ny ($(is_ny_valid ? "optimal" : "non-optimal "))")
         println("nz = $nz ($(is_nz_valid ? "optimal" : "non-optimal "))")
         println("total_ranks = $total_ranks ($(is_ranks_valid ? "optimal" : "non-optimal"))")
-        
+
         println("\nSuggested optimal configuration:")
         println("nx = $(get(suggestions, "nx", nx))")
         println("ny = $(get(suggestions, "ny", ny))")
         println("nz = $(get(suggestions, "nz", nz))")
         println("total_ranks = $(get(suggestions, "total_ranks", total_ranks))")
     end
-    
+
     return
 end
 
@@ -1709,20 +1713,20 @@ Sequence includes powers of 2 and their multiples by specified multipliers.
 function generate_valid_rank_sequence(max_value::Int = 196608, multipliers::Vector{Int} = [3])
     valid_ranks = Int[]
     power_of_2 = 2
-    
+
     while power_of_2 <= max_value
         push!(valid_ranks, power_of_2)
-        
+
         for m in multipliers
             multiple = m * (power_of_2 ÷ 2)
             if multiple <= max_value
                 push!(valid_ranks, multiple)
             end
         end
-        
+
         power_of_2 *= 2
     end
-    
+
     return sort(unique(valid_ranks))
 end
 
@@ -1753,7 +1757,7 @@ function update_nel_xyz_in_file(filename::String, nx::Int, ny::Int, nz::Int)
             lines[i] = "    nel_z   =  $nz"
         end
     end
-    open(filename, "w") do io
+    return open(filename, "w") do io
         for line in lines
             println(io, line)
         end
@@ -1796,21 +1800,21 @@ function insert_cpu_lines_after_nelz(filename::String, P)
     newlines = String[]
     inserted = false
     skip_next_lines = 0
-    
+
     for (i, line) in enumerate(lines)
         if skip_next_lines > 0
             skip_next_lines -= 1
             continue
         end
-        
+
         push!(newlines, line)
-        
+
         if !inserted && occursin("nel_z", line)
             # Check if CPU lines already exist in the next few lines
             if i + 4 <= length(lines) &&
-               occursin("cpu_x", lines[i+1]) &&
-               occursin("cpu_y", lines[i+2]) &&
-               occursin("cpu_z", lines[i+3])
+                    occursin("cpu_x", lines[i + 1]) &&
+                    occursin("cpu_y", lines[i + 2]) &&
+                    occursin("cpu_z", lines[i + 3])
                 # Replace existing lines
                 newlines[end] = ""  # Replace the current line with empty line
                 push!(newlines, "# Number of processes")
@@ -1829,10 +1833,10 @@ function insert_cpu_lines_after_nelz(filename::String, P)
             inserted = true
         end
     end
-    
-    open(filename, "w") do io
+
+    return open(filename, "w") do io
         for line in newlines
             println(io, line)
         end
     end
-end
\ No newline at end of file
+end
diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl
index 4b0baf2..ec934e8 100644
--- a/src/Setup_geometry.jl
+++ b/src/Setup_geometry.jl
@@ -323,32 +323,34 @@ end
     cell = false     )
 Add box function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing 
 """
-function  add_box!(Phase, Temp, Grid::AbstractGeneralGrid,           # required input
-                bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}};    # limits of the box
-                Origin = nothing, StrikeAngle = 0, DipAngle = 0,     # origin & dip/strike
-                phase = ConstantPhase(1),                            # Sets the phase number(s) in the box
-                T = nothing,                                         # Sets the thermal structure (various functions are available)
-                segments = nothing,                                  # Allows defining multiple ridge segments
-                cell = false                                         # if true, Phase and Temp are defined on cell centers
-                ) 
+function add_box!(
+        Phase, Temp, Grid::AbstractGeneralGrid,           # required input
+        bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}};    # limits of the box
+        Origin = nothing, StrikeAngle = 0, DipAngle = 0,     # origin & dip/strike
+        phase = ConstantPhase(1),                            # Sets the phase number(s) in the box
+        T = nothing,                                         # Sets the thermal structure (various functions are available)
+        segments = nothing,                                  # Allows defining multiple ridge segments
+        cell = false                                         # if true, Phase and Temp are defined on cell centers
+    )
 
     if isempty(bounds)
 
         return nothing
-        
+
     else
-        
-        xlim=(Tuple(bounds[1])) 
-        ylim=(Tuple(bounds[2])) 
-        zlim=(Tuple(bounds[3]))
-
-        add_box!( 
-                Phase, Temp, Grid;       # required input
-                xlim     = xlim, ylim = ylim, zlim = zlim,     # limits of the box
-                Origin   = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle,      # origin & dip/strike
-                phase    = phase,                          # Sets the phase number(s) in the box
-                T        = T,                                  # Sets the thermal structure (various functions are available)
-                cell     = cell )
+
+        xlim = (Tuple(bounds[1]))
+        ylim = (Tuple(bounds[2]))
+        zlim = (Tuple(bounds[3]))
+
+        add_box!(
+            Phase, Temp, Grid;       # required input
+            xlim = xlim, ylim = ylim, zlim = zlim,     # limits of the box
+            Origin = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle,      # origin & dip/strike
+            phase = phase,                          # Sets the phase number(s) in the box
+            T = T,                                  # Sets the thermal structure (various functions are available)
+            cell = cell
+        )
     end
 
 end
@@ -818,24 +820,26 @@ end
     cell = false     )
 Add polygon function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing 
 """
-function  add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid,           # required input
-                bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}};    # limits of the box
-                phase = ConstantPhase(1),                            # Sets the phase number(s) in the box
-                T = nothing,                                         # Sets the thermal structure (various functions are available)
-                cell = false                                         # if true, Phase and Temp are defined on cell centers
-                ) 
-
-    if !isempty(bounds)        
-        xlim=(Tuple(bounds[1])) 
-        ylim=(Tuple(bounds[2])) 
-        zlim=(Tuple(bounds[3]))
-
-        add_polygon!( 
-                Phase, Temp, Grid;       # required input
-                xlim     = xlim, ylim = ylim, zlim = zlim,     # limits of the box
-                phase    = phase,                              # Sets the phase number(s) in the box
-                T        = T,                                  # Sets the thermal structure (various functions are available)
-                cell     = cell )
+function add_polygon!(
+        Phase, Temp, Grid::AbstractGeneralGrid,           # required input
+        bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}};    # limits of the box
+        phase = ConstantPhase(1),                            # Sets the phase number(s) in the box
+        T = nothing,                                         # Sets the thermal structure (various functions are available)
+        cell = false                                         # if true, Phase and Temp are defined on cell centers
+    )
+
+    return if !isempty(bounds)
+        xlim = (Tuple(bounds[1]))
+        ylim = (Tuple(bounds[2]))
+        zlim = (Tuple(bounds[3]))
+
+        add_polygon!(
+            Phase, Temp, Grid;       # required input
+            xlim = xlim, ylim = ylim, zlim = zlim,     # limits of the box
+            phase = phase,                              # Sets the phase number(s) in the box
+            T = T,                                  # Sets the thermal structure (various functions are available)
+            cell = cell
+        )
     end
 
 end
@@ -936,27 +940,29 @@ end
     cell = false     )
 Add plate function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing 
 """
-function  add_plate!(Phase, Temp, Grid::AbstractGeneralGrid,           # required input
-                bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}};    # limits of the box
-                phase = ConstantPhase(1),                            # Sets the phase number(s) in the box
-                T = nothing,                                         # Sets the thermal structure (various functions are available)
-                segments = nothing,                                  # Allows defining multiple ridge segments
-                cell = false                                         # if true, Phase and Temp are defined on cell centers
-                ) 
-
-    if !isempty(bounds)        
-        xlim=(Tuple(bounds[1])) 
-        ylim=(Tuple(bounds[2])) 
-        zlim=(Tuple(bounds[3]))
-
-        add_plate!( 
-                Phase, Temp, Grid;       # required input
-                xlim     = xlim, ylim = ylim, zlim = zlim,     # limits of the box
-                Origin   = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle,      # origin & dip/strike
-                phase    = phase,                          # Sets the phase number(s) in the box
-                T        = T,                                  # Sets the thermal structure (various functions are available)
-                segments = segments,                     # Allows defining multiple ridge segments
-                cell     = cell )
+function add_plate!(
+        Phase, Temp, Grid::AbstractGeneralGrid,           # required input
+        bounds::Union{Vector{Any}, AbstractVector{<:AbstractVector{<:Real}}};    # limits of the box
+        phase = ConstantPhase(1),                            # Sets the phase number(s) in the box
+        T = nothing,                                         # Sets the thermal structure (various functions are available)
+        segments = nothing,                                  # Allows defining multiple ridge segments
+        cell = false                                         # if true, Phase and Temp are defined on cell centers
+    )
+
+    return if !isempty(bounds)
+        xlim = (Tuple(bounds[1]))
+        ylim = (Tuple(bounds[2]))
+        zlim = (Tuple(bounds[3]))
+
+        add_plate!(
+            Phase, Temp, Grid;       # required input
+            xlim = xlim, ylim = ylim, zlim = zlim,     # limits of the box
+            Origin = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle,      # origin & dip/strike
+            phase = phase,                          # Sets the phase number(s) in the box
+            T = T,                                  # Sets the thermal structure (various functions are available)
+            segments = segments,                     # Allows defining multiple ridge segments
+            cell = cell
+        )
     end
 end
 
diff --git a/test/test_lamem.jl b/test/test_lamem.jl
index a120925..189effa 100644
--- a/test/test_lamem.jl
+++ b/test/test_lamem.jl
@@ -47,73 +47,73 @@ save_LaMEM_markers_parallel(Model3D, verbose = false)
 save_LaMEM_markers_parallel(Model3D, PartitioningFile = PartitioningFile, verbose = false)
 
 # Get partitioning without generating partitioning file
-n_ranks = 8; nx=128; ny=64; nz=32;
-xcoords_ans=[-35.35034129014803, -19.867446943762857, -4.3845525973776835, 11.098341749007488, 26.58123609539267];
-ycoords_ans=[-24.510578171895816, 5.122856149231547, 34.75629047035891];
-zcoords_ans=[-6.4, 6.4];
+n_ranks = 8; nx = 128; ny = 64; nz = 32;
+xcoords_ans = [-35.35034129014803, -19.867446943762857, -4.3845525973776835, 11.098341749007488, 26.58123609539267];
+ycoords_ans = [-24.510578171895816, 5.122856149231547, 34.75629047035891];
+zcoords_ans = [-6.4, 6.4];
 nProcX_ans = 4; nProcY_ans = 2; nProcZ_ans = 1;
 nNodeX_ans = 129; nNodeY_ans = 65; nNodeZ_ans = 33;
-P         = setup_model_domain([extrema(xcoords_ans)...], [extrema(ycoords_ans)...], [extrema(zcoords_ans)...], nx, ny, nz, n_ranks)
-@test isapprox(nProcX_ans, P.nProcX; atol=1.0e-8)
-@test isapprox(nProcY_ans, P.nProcY; atol=1.0e-8)
-@test isapprox(nProcZ_ans, P.nProcZ; atol=1.0e-8)
-@test isapprox(xcoords_ans, P.xc; atol=1.0e-8)
-@test isapprox(ycoords_ans, P.yc; atol=1.0e-8)
-@test isapprox(zcoords_ans, P.zc; atol=1.0e-8)
-@test isapprox(nNodeX_ans, P.nNodeX; atol=1.0e-8)
-@test isapprox(nNodeY_ans, P.nNodeY; atol=1.0e-8)
-@test isapprox(nNodeZ_ans, P.nNodeZ; atol=1.0e-8)
-
-n_ranks = 128; nx=256; ny=1; nz=128;
-xcoords_ans=[-35.35034129014803, -31.479617703551735, -27.608894116955444, -23.73817053035915, -19.867446943762857, -15.996723357166562, -12.125999770570267, -8.255276183973974, -4.384552597377681, -0.5138290107813872, 3.3568945758149065, 7.227618162411201, 11.098341749007494, 14.969065335603787, 18.839788922200082, 22.710512508796374, 26.58123609539267];
-ycoords_ans=[-0.1, 0.1];
-zcoords_ans=[-6.4, -4.8, -3.2, -1.6, 0.0, 1.6, 3.2, 4.8, 6.4];
+P = setup_model_domain([extrema(xcoords_ans)...], [extrema(ycoords_ans)...], [extrema(zcoords_ans)...], nx, ny, nz, n_ranks)
+@test isapprox(nProcX_ans, P.nProcX; atol = 1.0e-8)
+@test isapprox(nProcY_ans, P.nProcY; atol = 1.0e-8)
+@test isapprox(nProcZ_ans, P.nProcZ; atol = 1.0e-8)
+@test isapprox(xcoords_ans, P.xc; atol = 1.0e-8)
+@test isapprox(ycoords_ans, P.yc; atol = 1.0e-8)
+@test isapprox(zcoords_ans, P.zc; atol = 1.0e-8)
+@test isapprox(nNodeX_ans, P.nNodeX; atol = 1.0e-8)
+@test isapprox(nNodeY_ans, P.nNodeY; atol = 1.0e-8)
+@test isapprox(nNodeZ_ans, P.nNodeZ; atol = 1.0e-8)
+
+n_ranks = 128; nx = 256; ny = 1; nz = 128;
+xcoords_ans = [-35.35034129014803, -31.479617703551735, -27.608894116955444, -23.73817053035915, -19.867446943762857, -15.996723357166562, -12.125999770570267, -8.255276183973974, -4.384552597377681, -0.5138290107813872, 3.3568945758149065, 7.227618162411201, 11.098341749007494, 14.969065335603787, 18.839788922200082, 22.710512508796374, 26.58123609539267];
+ycoords_ans = [-0.1, 0.1];
+zcoords_ans = [-6.4, -4.8, -3.2, -1.6, 0.0, 1.6, 3.2, 4.8, 6.4];
 nProcX_ans = 16; nProcY_ans = 1; nProcZ_ans = 8;
 nNodeX_ans = 257; nNodeY_ans = 2; nNodeZ_ans = 129;
-P         = setup_model_domain([extrema(xcoords_ans)...], [extrema(ycoords_ans)...], [extrema(zcoords_ans)...], nx, ny, nz, n_ranks)
-@test isapprox(nProcX_ans, P.nProcX; atol=1.0e-8)
-@test isapprox(nProcY_ans, P.nProcY; atol=1.0e-8)
-@test isapprox(nProcZ_ans, P.nProcZ; atol=1.0e-8)
-@test isapprox(xcoords_ans, P.xc; atol=1.0e-8)
-@test isapprox(ycoords_ans, P.yc; atol=1.0e-8)
-@test isapprox(zcoords_ans, P.zc; atol=1.0e-8)
-@test isapprox(nNodeX_ans, P.nNodeX; atol=1.0e-8)
-@test isapprox(nNodeY_ans, P.nNodeY; atol=1.0e-8)
-@test isapprox(nNodeZ_ans, P.nNodeZ; atol=1.0e-8)
-
-n_ranks = 2048; nx=512; ny=2048; nz=512;
-xcoords_ans=[-35.35034129014803, -27.60889411695544, -19.867446943762857, -12.12599977057027, -4.3845525973776835, 3.356894575814903, 11.098341749007488, 18.839788922200068, 26.58123609539267 ];
-ycoords_ans=[-24.510578171895816, -22.658488526825355, -20.806398881754898, -18.954309236684434, -17.102219591613977, -15.250129946543517, -13.398040301473054, -11.545950656402596,  -9.693861011332135,  -7.841771366261674];
-zcoords_ans=[-6.4, -4.800000000000001, -3.2, -1.6, 0.0, 1.6, 3.2, 4.800000000000001, 6.4];
+P = setup_model_domain([extrema(xcoords_ans)...], [extrema(ycoords_ans)...], [extrema(zcoords_ans)...], nx, ny, nz, n_ranks)
+@test isapprox(nProcX_ans, P.nProcX; atol = 1.0e-8)
+@test isapprox(nProcY_ans, P.nProcY; atol = 1.0e-8)
+@test isapprox(nProcZ_ans, P.nProcZ; atol = 1.0e-8)
+@test isapprox(xcoords_ans, P.xc; atol = 1.0e-8)
+@test isapprox(ycoords_ans, P.yc; atol = 1.0e-8)
+@test isapprox(zcoords_ans, P.zc; atol = 1.0e-8)
+@test isapprox(nNodeX_ans, P.nNodeX; atol = 1.0e-8)
+@test isapprox(nNodeY_ans, P.nNodeY; atol = 1.0e-8)
+@test isapprox(nNodeZ_ans, P.nNodeZ; atol = 1.0e-8)
+
+n_ranks = 2048; nx = 512; ny = 2048; nz = 512;
+xcoords_ans = [-35.35034129014803, -27.60889411695544, -19.867446943762857, -12.12599977057027, -4.3845525973776835, 3.356894575814903, 11.098341749007488, 18.839788922200068, 26.58123609539267];
+ycoords_ans = [-24.510578171895816, -22.658488526825355, -20.806398881754898, -18.954309236684434, -17.102219591613977, -15.250129946543517, -13.398040301473054, -11.545950656402596, -9.693861011332135, -7.841771366261674];
+zcoords_ans = [-6.4, -4.800000000000001, -3.2, -1.6, 0.0, 1.6, 3.2, 4.800000000000001, 6.4];
 nProcX_ans = 8; nProcY_ans = 32; nProcZ_ans = 8;
 nNodeX_ans = 513; nNodeY_ans = 2049; nNodeZ_ans = 513;
-P         = setup_model_domain([-35.35034129014803,26.58123609539267], [-24.510578171895816,34.75629047035891], [-6.4,6.4], nx, ny, nz, n_ranks)
-@test isapprox(nProcX_ans, P.nProcX; atol=1.0e-8)
-@test isapprox(nProcY_ans, P.nProcY; atol=1.0e-8)
-@test isapprox(nProcZ_ans, P.nProcZ; atol=1.0e-8)
-@test isapprox(xcoords_ans, P.xc; atol=1.0e-8)
-@test isapprox(ycoords_ans, P.yc[1:10]; atol=1.0e-8)
-@test isapprox(zcoords_ans, P.zc; atol=1.0e-8)
-@test isapprox(nNodeX_ans, P.nNodeX; atol=1.0e-8)
-@test isapprox(nNodeY_ans, P.nNodeY; atol=1.0e-8)
-@test isapprox(nNodeZ_ans, P.nNodeZ; atol=1.0e-8)
-
-n_ranks = 32768; nx=2048; ny=512; nz=1024;
-xcoords_ans=[-35.35034129014803, -34.38266039349895, -33.41497949684988, -32.4472986002008, -31.47961770355173, -30.511936806902664, -29.54425591025359, -28.576575013604515, -27.60889411695544, -26.64121322030637 ];
-ycoords_ans=[ -24.510578171895816, -20.806398881754898, -17.102219591613977, -13.398040301473054,  -9.693861011332135,  -5.989681721191215,  -2.285502431050293,   1.418676859090624,   5.122856149231547,   8.82703543937247];
-zcoords_ans=[-6.4, -6.0, -5.6000000000000005, -5.2, -4.800000000000001, -4.4, -4.0, -3.6, -3.2, -2.8000000000000003];
+P = setup_model_domain([-35.35034129014803, 26.58123609539267], [-24.510578171895816, 34.75629047035891], [-6.4, 6.4], nx, ny, nz, n_ranks)
+@test isapprox(nProcX_ans, P.nProcX; atol = 1.0e-8)
+@test isapprox(nProcY_ans, P.nProcY; atol = 1.0e-8)
+@test isapprox(nProcZ_ans, P.nProcZ; atol = 1.0e-8)
+@test isapprox(xcoords_ans, P.xc; atol = 1.0e-8)
+@test isapprox(ycoords_ans, P.yc[1:10]; atol = 1.0e-8)
+@test isapprox(zcoords_ans, P.zc; atol = 1.0e-8)
+@test isapprox(nNodeX_ans, P.nNodeX; atol = 1.0e-8)
+@test isapprox(nNodeY_ans, P.nNodeY; atol = 1.0e-8)
+@test isapprox(nNodeZ_ans, P.nNodeZ; atol = 1.0e-8)
+
+n_ranks = 32768; nx = 2048; ny = 512; nz = 1024;
+xcoords_ans = [-35.35034129014803, -34.38266039349895, -33.41497949684988, -32.4472986002008, -31.47961770355173, -30.511936806902664, -29.54425591025359, -28.576575013604515, -27.60889411695544, -26.64121322030637];
+ycoords_ans = [-24.510578171895816, -20.806398881754898, -17.102219591613977, -13.398040301473054, -9.693861011332135, -5.989681721191215, -2.285502431050293, 1.418676859090624, 5.122856149231547, 8.82703543937247];
+zcoords_ans = [-6.4, -6.0, -5.6000000000000005, -5.2, -4.800000000000001, -4.4, -4.0, -3.6, -3.2, -2.8000000000000003];
 nProcX_ans = 64; nProcY_ans = 16; nProcZ_ans = 32;
 nNodeX_ans = 2049; nNodeY_ans = 513; nNodeZ_ans = 1025;
-P         = setup_model_domain([-35.35034129014803,26.58123609539267], [-24.510578171895816,34.75629047035891], [-6.4,6.4], nx, ny, nz, n_ranks)
-@test isapprox(nProcX_ans, P.nProcX; atol=1.0e-8)
-@test isapprox(nProcY_ans, P.nProcY; atol=1.0e-8)
-@test isapprox(nProcZ_ans, P.nProcZ; atol=1.0e-8)
-@test isapprox(xcoords_ans, P.xc[1:10]; atol=1.0e-8)
-@test isapprox(ycoords_ans, P.yc[1:10]; atol=1.0e-8)
-@test isapprox(zcoords_ans, P.zc[1:10]; atol=1.0e-8)
-@test isapprox(nNodeX_ans, P.nNodeX; atol=1.0e-8)
-@test isapprox(nNodeY_ans, P.nNodeY; atol=1.0e-8)
-@test isapprox(nNodeZ_ans, P.nNodeZ; atol=1.0e-8)
+P = setup_model_domain([-35.35034129014803, 26.58123609539267], [-24.510578171895816, 34.75629047035891], [-6.4, 6.4], nx, ny, nz, n_ranks)
+@test isapprox(nProcX_ans, P.nProcX; atol = 1.0e-8)
+@test isapprox(nProcY_ans, P.nProcY; atol = 1.0e-8)
+@test isapprox(nProcZ_ans, P.nProcZ; atol = 1.0e-8)
+@test isapprox(xcoords_ans, P.xc[1:10]; atol = 1.0e-8)
+@test isapprox(ycoords_ans, P.yc[1:10]; atol = 1.0e-8)
+@test isapprox(zcoords_ans, P.zc[1:10]; atol = 1.0e-8)
+@test isapprox(nNodeX_ans, P.nNodeX; atol = 1.0e-8)
+@test isapprox(nNodeY_ans, P.nNodeY; atol = 1.0e-8)
+@test isapprox(nNodeZ_ans, P.nNodeZ; atol = 1.0e-8)
 
 # Test creating model setups
 Grid = read_LaMEM_inputfile("test_files/Subduction2D_FreeSlip_Particles_Linear_DirectSolver.dat")

[Copilot]

Pull Request Overview

This PR updates routines for LaMEM by enhancing geometry function interfaces to accept bounds as vectors and expanding test coverage.

• Adds tests in test_setup_geometry.jl for add_box! handling empty and non-empty bounds
• Enhances test_lamem.jl validations for setup_model_domain functions and grid partitioning
• Introduces new bounds-based functions in src/Setup_geometry.jl for box, polygon, and plate addition

Reviewed Changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated no comments.

File	Description
test/test_setup_geometry.jl	Added tests for add_box! with both empty and filled bounds to verify behavior
test/test_lamem.jl	Introduces additional tests checking domain partitioning and grid coordinate computations
src/Setup_geometry.jl	Adds new functions for bounds-based addition of boxes, polygons, and plates; potential bug in add_plate!
.vscode/settings.json	Contains a trivial empty settings file

Comments suppressed due to low confidence (2)

test/test_setup_geometry.jl:17

• The variable 'Data' appears to be undefined. Consider using a defined variable, such as 'Grid', if that was intended.

Phases = zeros(Int32, size(Data));

src/Setup_geometry.jl:965

• The variables 'Origin', 'StrikeAngle', and 'DipAngle' are not declared in the function parameters for add_plate!. Either remove these references or add them to the function signature if needed.

Origin   = Origin, StrikeAngle = StrikeAngle, DipAngle = DipAngle,

[boriskaus]

please remove the accidentaly added file; afterwards it seems fine to be merged

[albert-de-montserrat]

Could please format the code as well? thanks!