From bb0e7a457a60d343cfcd054b0e451d98f1904725 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 8 Oct 2024 17:49:54 +0200 Subject: [PATCH 01/36] store 1D vectors instead of ranges --- src/data_types.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/data_types.jl b/src/data_types.jl index 9eddd458..b3786d2d 100644 --- a/src/data_types.jl +++ b/src/data_types.jl @@ -1095,8 +1095,8 @@ struct CartGrid{FT, D} <: AbstractGeneralGrid L :: NTuple{D,FT} # Domain size min :: NTuple{D,FT} # start of the grid in every direction max :: NTuple{D,FT} # end of the grid in every direction - coord1D :: NTuple{D,StepRangeLen{FT}} # Tuple with 1D vectors in all directions - coord1D_cen :: NTuple{D,StepRangeLen{FT}} # Tuple with 1D vectors of center points in all directions + coord1D :: NTuple{D,Vector{FT}} # Tuple with 1D vectors in all directions + coord1D_cen :: NTuple{D,Vector{FT}} # Tuple with 1D vectors of center points in all directions end size(d::CartGrid) = d.N @@ -1202,13 +1202,13 @@ function create_CartGrid(; # Generate 1D coordinate arrays of vertices in all directions coord1D=() for idim=1:dim - coord1D = (coord1D..., range(X₁[idim], Xₙ[idim]; length = N[idim] )) + coord1D = (coord1D..., Vector(range(X₁[idim], Xₙ[idim]; length = N[idim] ))) end # Generate 1D coordinate arrays centers in all directionbs coord1D_cen=() for idim=1:dim - coord1D_cen = (coord1D_cen..., range(X₁[idim]+Δ[idim]/2, Xₙ[idim]-Δ[idim]/2; length = N[idim]-1 )) + coord1D_cen = (coord1D_cen..., Vector(range(X₁[idim]+Δ[idim]/2, Xₙ[idim]-Δ[idim]/2; length = N[idim]-1 ))) end ConstantΔ = true; From 53275a725ed70cb3530da479a0ce390ccc44f03e Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 8 Oct 2024 19:11:33 +0200 Subject: [PATCH 02/36] add Chmy interface and multiple dispatch for most functions --- Project.toml | 71 ++++++++++++++++++++++++----------------------- ext/Chmy_utils.jl | 69 +++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 4 ++- test/test_Chmy.jl | 40 ++++++++++++++++++++++++++ 4 files changed, 149 insertions(+), 35 deletions(-) create mode 100644 ext/Chmy_utils.jl create mode 100644 test/test_Chmy.jl diff --git a/Project.toml b/Project.toml index 53d380c7..dc986ec9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,34 @@ +authors = ["Boris Kaus", "Marcel Thielmann"] name = "GeophysicalModelGenerator" uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742" -authors = ["Boris Kaus", "Marcel Thielmann"] -version = "0.7.7" +version = "0.7.8" + +[compat] +Chmy = "0.1.20" +Colors = "0.9 - 0.12" +Downloads = "1" +FFMPEG = "0.4" +FileIO = "1" +GDAL_jll = "300.900.0 - 301.901.0" +GLMakie = "0.8, 0.9, 0.10" +GMT = "1.0 - 1.14" +GeoParams = "0.2 - 0.6" +Geodesy = "1" +GeometryBasics = "0.1 - 0.4" +Glob = "1.2 - 1.3" +GridapGmsh = "0.5 - 0.7" +ImageIO = "0.1 - 0.6" +Interpolations = "0.14, 0.15" +JLD2 = "0.4, 0.5" +LightXML = "0.8, 0.9" +MeshIO = "0.1 - 0.4" +NearestNeighbors = "0.2 - 0.4" +Parameters = "0.9 - 0.12" +SpecialFunctions = "1.0, 2" +StaticArrays = "1" +WhereTheWaterFlows = "0.10, 0.11" +WriteVTK = "1" +julia = "1.9" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" @@ -31,41 +58,11 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" WhereTheWaterFlows = "ea906314-1493-4d22-a0af-f886a20c9fba" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" -[weakdeps] -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" -GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" -GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8" - [extensions] GLMakie_Visualisation = "GLMakie" GMT_utils = "GMT" Gmsh_utils = "GridapGmsh" - -[compat] -Colors = "0.9 - 0.12" -Downloads = "1" -FFMPEG = "0.4" -FileIO = "1" -GDAL_jll = "300.900.0 - 301.900.0" -GLMakie = "0.8, 0.9, 0.10" -GMT = "1.0 - 1.14" -GeoParams = "0.2 - 0.6" -Geodesy = "1" -GeometryBasics = "0.1 - 0.4" -Glob = "1.2 - 1.3" -GridapGmsh = "0.5 - 0.7" -ImageIO = "0.1 - 0.6" -Interpolations = "0.14, 0.15" -JLD2 = "0.4, 0.5" -LightXML = "0.8, 0.9" -MeshIO = "0.1 - 0.4" -NearestNeighbors = "0.2 - 0.4" -Parameters = "0.9 - 0.12" -SpecialFunctions = "1.0, 2" -StaticArrays = "1" -WhereTheWaterFlows = "0.10" -WriteVTK = "1" -julia = "1.9" +Chmy_utils = "Chmy" [extras] GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" @@ -74,4 +71,10 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "GMT", "StableRNGs", "GridapGmsh"] +test = ["Test", "GMT", "StableRNGs", "GridapGmsh", "Chmy"] + +[weakdeps] +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" +GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8" +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" diff --git a/ext/Chmy_utils.jl b/ext/Chmy_utils.jl new file mode 100644 index 00000000..4faca7ad --- /dev/null +++ b/ext/Chmy_utils.jl @@ -0,0 +1,69 @@ +#helper functions to make GMG work with Chmy grids and fields +module Chmy_utils + +using Chmy, Chmy.Grids, Chmy.Fields + +import GeophysicalModelGenerator: create_CartGrid, CartGrid +import GeophysicalModelGenerator: add_box!, add_sphere!, add_ellipsoid!, add_cylinder! +import GeophysicalModelGenerator: add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano! + +println("Loading Chmy-GMG tools") + +""" + CartGrid = create_CartGrid(grid::StructuredGrid) + +Creates a GMG `CartGrid` data structure from a `Chmy` grid object +""" +function create_CartGrid(grid::StructuredGrid) + + coord1D = Vector.(coords(grid, Vertex())) + coord1D_cen = Vector.(coords(grid, Center())) + N = length.(coord1D) + L = extent(grid, Vertex()) + X₁ = origin(grid, Vertex()) + Xₙ = X₁ .+ L + Δ = spacing(grid) + ConstantΔ = false; + if isa(typeof(grid), UniformGrid) + ConstantΔ = true + end + + return CartGrid(ConstantΔ,N,Δ,L,X₁,Xₙ,coord1D, coord1D_cen) +end + +# all functions to be ported +function_names = (:add_box!, :add_sphere!, :add_ellipsoid!, :add_cylinder!, :add_layer!, :add_polygon!, :add_slab!, :add_stripes!, :add_volcano!) + +for fn in function_names + + @eval begin + """ + $($fn)( Phase::Field, + Temp::Field, + Grid::StructuredGrid; # required input + kwargs...) + + Sets `$($fn)` function for `Chmy` fields and grids. + """ + function $fn( Phase::Field, + Temp::Field, + Grid::StructuredGrid; # required input + kwargs...) + + CartGrid = create_CartGrid(Grid) + + cell = false + if all(location(Phase).==Center()) + cell = true + end + + return ($fn)(Phase, Temp, CartGrid; cell=cell, kwargs...) + end + end + + + +end + + +end # end of module diff --git a/test/runtests.jl b/test/runtests.jl index e6438262..d9c61b84 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -80,7 +80,9 @@ using Test @testset "ASAGI_IO" begin include("test_ASAGI_IO.jl") end - + @testset "Chmy" begin + include("test_Chmy.jl") + end end # Cleanup diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl new file mode 100644 index 00000000..e3bc6e28 --- /dev/null +++ b/test/test_Chmy.jl @@ -0,0 +1,40 @@ +using Test, GeophysicalModelGenerator + +using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.Fields +using KernelAbstractions # for backend-agnostic kernels + +backend = CPU() +arch = Arch(backend) + +# 3D test +lx, ly, lz = 10.0, 11.0, 12.0 +nx, ny, nz = 10,11,12 +grid = UniformGrid(arch; + origin=(-lx/2, -ly/2, -lz/2), + extent=(lx, ly, lz), + dims=(nx, ny, nz)) + +# create field +Temp_C = Field(backend, grid, Center(), Float64; halo=1) +Phases_C = Field(backend, grid, Center(), Int32; halo=1) +Temp_V = Field(backend, grid, Vertex(), Float64; halo=1) +Phases_V = Field(backend, grid, Vertex(), Int32; halo=1) + +# grid: +CartGrid = create_CartGrid(grid) + +@test sum.(CartGrid.coord1D) == (0.0, 0.0, 0.0) + +# test add_box! directly. Note that this requires you to specify a "cell" keyword for Center() localtions +add_box!(Phases_C,Temp_C,CartGrid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(3), cell=true) +@test extrema(Phases_C) == (0,3) + +add_box!(Phases_V,Temp_V,CartGrid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(3)) +@test extrema(Phases_V) == (0,3) + +# multiple dispatch functions +add_box!(Phases_C,Temp_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2)) +@test extrema(Phases_C) == (0,2) + +add_box!(Phases_V,Temp_V,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2)) +@test extrema(Phases_V) == (0,2) From 884272fdb1e17fd501adebe4f63deada19ab968f Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 8 Oct 2024 23:17:30 +0200 Subject: [PATCH 03/36] missing dependencies --- Project.toml | 74 +++++++++++++++++++++++++++------------------------- 1 file changed, 38 insertions(+), 36 deletions(-) diff --git a/Project.toml b/Project.toml index dc986ec9..1517f95b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,35 +1,8 @@ -authors = ["Boris Kaus", "Marcel Thielmann"] name = "GeophysicalModelGenerator" uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742" +authors = ["Boris Kaus", "Marcel Thielmann"] version = "0.7.8" -[compat] -Chmy = "0.1.20" -Colors = "0.9 - 0.12" -Downloads = "1" -FFMPEG = "0.4" -FileIO = "1" -GDAL_jll = "300.900.0 - 301.901.0" -GLMakie = "0.8, 0.9, 0.10" -GMT = "1.0 - 1.14" -GeoParams = "0.2 - 0.6" -Geodesy = "1" -GeometryBasics = "0.1 - 0.4" -Glob = "1.2 - 1.3" -GridapGmsh = "0.5 - 0.7" -ImageIO = "0.1 - 0.6" -Interpolations = "0.14, 0.15" -JLD2 = "0.4, 0.5" -LightXML = "0.8, 0.9" -MeshIO = "0.1 - 0.4" -NearestNeighbors = "0.2 - 0.4" -Parameters = "0.9 - 0.12" -SpecialFunctions = "1.0, 2" -StaticArrays = "1" -WhereTheWaterFlows = "0.10, 0.11" -WriteVTK = "1" -julia = "1.9" - [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" @@ -58,23 +31,52 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" WhereTheWaterFlows = "ea906314-1493-4d22-a0af-f886a20c9fba" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" +[weakdeps] +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" +GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8" + [extensions] +Chmy_utils = "Chmy" GLMakie_Visualisation = "GLMakie" GMT_utils = "GMT" Gmsh_utils = "GridapGmsh" -Chmy_utils = "Chmy" + +[compat] +Chmy = "0.1.20" +Colors = "0.9 - 0.12" +Downloads = "1" +FFMPEG = "0.4" +FileIO = "1" +GDAL_jll = "300.900.0 - 301.901.0" +GLMakie = "0.8, 0.9, 0.10" +GMT = "1.0 - 1.14" +GeoParams = "0.2 - 0.6" +Geodesy = "1" +GeometryBasics = "0.1 - 0.4" +Glob = "1.2 - 1.3" +GridapGmsh = "0.5 - 0.7" +ImageIO = "0.1 - 0.6" +Interpolations = "0.14, 0.15" +JLD2 = "0.4, 0.5" +LightXML = "0.8, 0.9" +MeshIO = "0.1 - 0.4" +NearestNeighbors = "0.2 - 0.4" +Parameters = "0.9 - 0.12" +SpecialFunctions = "1.0, 2" +StaticArrays = "1" +WhereTheWaterFlows = "0.10, 0.11" +WriteVTK = "1" +julia = "1.9" [extras] +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "GMT", "StableRNGs", "GridapGmsh", "Chmy"] - -[weakdeps] -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" -GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54" -GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8" -Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" +test = ["Test", "GMT", "StableRNGs", "GridapGmsh", "Chmy","KernelAbstractions"] From 4abc21fbabaf6f36f1bc8799625355d251bab14e Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Tue, 8 Oct 2024 23:22:46 +0200 Subject: [PATCH 04/36] spellcheck --- test/test_Chmy.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl index e3bc6e28..a7f7a747 100644 --- a/test/test_Chmy.jl +++ b/test/test_Chmy.jl @@ -25,7 +25,7 @@ CartGrid = create_CartGrid(grid) @test sum.(CartGrid.coord1D) == (0.0, 0.0, 0.0) -# test add_box! directly. Note that this requires you to specify a "cell" keyword for Center() localtions +# test add_box! directly. Note that this requires you to specify a "cell" keyword for Center() locations add_box!(Phases_C,Temp_C,CartGrid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(3), cell=true) @test extrema(Phases_C) == (0,3) From 0cf7795696b2e09d3cee56863286dea01eba85ac Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 09:36:33 +0200 Subject: [PATCH 05/36] add functionality for above/belowSurface --- ext/Chmy_utils.jl | 33 +++++++++++++++++++++++++++++++-- src/surface_functions.jl | 12 ++++++++---- test/gmt.history | 2 +- test/test_Chmy.jl | 16 ++++++++++++++++ 4 files changed, 56 insertions(+), 7 deletions(-) diff --git a/ext/Chmy_utils.jl b/ext/Chmy_utils.jl index 4faca7ad..675d2bf5 100644 --- a/ext/Chmy_utils.jl +++ b/ext/Chmy_utils.jl @@ -3,10 +3,11 @@ module Chmy_utils using Chmy, Chmy.Grids, Chmy.Fields -import GeophysicalModelGenerator: create_CartGrid, CartGrid +import GeophysicalModelGenerator: create_CartGrid, CartGrid, CartData import GeophysicalModelGenerator: add_box!, add_sphere!, add_ellipsoid!, add_cylinder! import GeophysicalModelGenerator: add_layer!, add_polygon!, add_slab!, add_stripes!, add_volcano! - +import GeophysicalModelGenerator: above_surface, below_surface + println("Loading Chmy-GMG tools") """ @@ -61,8 +62,36 @@ for fn in function_names end end +end + + +# all functions to be ported +function_names = (:above_surface, :below_surface) + +for fn in function_names + + @eval begin + """ + $($fn)( Grid::StructuredGrid, field::Field, DataSurface_Cart::CartData; kwargs...) + Sets `$($fn)` function for `Chmy` grids and the field `field` whioch can be either on vertexes or centers + """ + function $fn( Grid::StructuredGrid, + field::Field, + DataSurface_Cart::CartData; + kwargs...) + + CartGrid = create_CartGrid(Grid) + + cell = false + if all(location(field).==Center()) + cell = true + end + return ($fn)(CartGrid, DataSurface_Cart; cell=cell, kwargs...) + end + end + end diff --git a/src/surface_functions.jl b/src/surface_functions.jl index 79cde200..17265abc 100644 --- a/src/surface_functions.jl +++ b/src/surface_functions.jl @@ -287,14 +287,18 @@ function above_surface(Data_Cart::Union{Q1Data,CartData}, DataSurface_Cart::Cart end """ - Above = above_surface(Grid::CartGrid, DataSurface_Cart::CartData; above=true) + Above = above_surface(Grid::CartGrid, DataSurface_Cart::CartData; above=true, cell=false) Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart` """ -function above_surface(Grid::CartGrid, DataSurface_Cart::CartData; above=true) +function above_surface(Grid::CartGrid, DataSurface_Cart::CartData; above=true, cell=false) - X,Y,Z = xyz_grid(Grid.coord1D...) - Data = CartData(Grid,(Z=Z,)) + if cell + X,Y,Z = xyz_grid(Grid.coord1D_cen...) + else + X,Y,Z = xyz_grid(Grid.coord1D...) + end + Data = CartData(X,Y,Z,(Z=Z,)) return above_surface(Data, DataSurface_Cart; above=above) end diff --git a/test/gmt.history b/test/gmt.history index 15bf1687..ebcde871 100644 --- a/test/gmt.history +++ b/test/gmt.history @@ -1,4 +1,4 @@ # GMT 6 Session common arguments shelf BEGIN GMT 6.6.0 -R 50/51/30/31 +R 6.5/7.3/50.2/50.6 END diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl index a7f7a747..fe7f7b28 100644 --- a/test/test_Chmy.jl +++ b/test/test_Chmy.jl @@ -38,3 +38,19 @@ add_box!(Phases_C,Temp_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2) add_box!(Phases_V,Temp_V,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2)) @test extrema(Phases_V) == (0,2) + +add_sphere!(Phases_V,Temp_V,grid, cen=(0,0,-1), radius=2.5, phase=ConstantPhase(3), T=ConstantTemp(800)) +@test extrema(Phases_V) == (0,3) +@test extrema(Temp_V) == (0.0,800.0) + +# test above/below surface intersection +Topo_cart = CartData(xyz_grid(-6:.2:6,-12:.2:13,0)); +ind = above_surface(grid, Phases_V, Topo_cart); +Phases_V[ind] .= 4; +@test extrema(Phases_V) == (0,3) + +ind = above_surface(grid, Phases_C, Topo_cart); +Phases_C[ind] .= 4; +@test extrema(Phases_V) == (0,4) + + From 5a3b2cdbaa8b3576fe8d32b25bad9de105dfeb18 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 09:43:43 +0200 Subject: [PATCH 06/36] spell check --- ext/Chmy_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/Chmy_utils.jl b/ext/Chmy_utils.jl index 675d2bf5..68bab739 100644 --- a/ext/Chmy_utils.jl +++ b/ext/Chmy_utils.jl @@ -74,7 +74,7 @@ for fn in function_names """ $($fn)( Grid::StructuredGrid, field::Field, DataSurface_Cart::CartData; kwargs...) - Sets `$($fn)` function for `Chmy` grids and the field `field` whioch can be either on vertexes or centers + Sets `$($fn)` function for `Chmy` grids and the field `field` which can be either on vertices or centers """ function $fn( Grid::StructuredGrid, field::Field, From 94daeeafe64bf2f381f549aac4fba425e3a28378 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 10:00:43 +0200 Subject: [PATCH 07/36] correct test --- test/test_Chmy.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl index fe7f7b28..901dee3d 100644 --- a/test/test_Chmy.jl +++ b/test/test_Chmy.jl @@ -47,10 +47,25 @@ add_sphere!(Phases_V,Temp_V,grid, cen=(0,0,-1), radius=2.5, phase=ConstantPhase Topo_cart = CartData(xyz_grid(-6:.2:6,-12:.2:13,0)); ind = above_surface(grid, Phases_V, Topo_cart); Phases_V[ind] .= 4; -@test extrema(Phases_V) == (0,3) +@test extrema(Phases_V) == (0,4) ind = above_surface(grid, Phases_C, Topo_cart); Phases_C[ind] .= 4; @test extrema(Phases_V) == (0,4) + + +# 2D test +lx, lz = 10.0, 12.0 +nx, nz = 10, 12 +grid = UniformGrid(arch; + origin=(-lx/2, -lz/2), + extent=(lx, lz), + dims=(nx, nz)) + +# create field +Temp2D_C = Field(backend, grid, Center()) +Phases2D_C = Field(backend, grid, Center(), Int32) + +add_box!(Phases2D_C,Temp2D_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2)) From 7a65cf23cbc7f6af59af95eda33a2f52e4aba94f Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 10:02:46 +0200 Subject: [PATCH 08/36] albert's suggestion --- src/data_types.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/data_types.jl b/src/data_types.jl index b3786d2d..fd7522f5 100644 --- a/src/data_types.jl +++ b/src/data_types.jl @@ -1089,14 +1089,14 @@ end Structure that holds data for an orthogonal cartesian grid, which can be described with 1D vectors """ struct CartGrid{FT, D} <: AbstractGeneralGrid - ConstantΔ :: Bool # Constant spacing (true in all cases for now) - N :: NTuple{D,Int} # Number of grid points in every direction - Δ :: NTuple{D,FT} # (constant) spacing in every direction - L :: NTuple{D,FT} # Domain size - min :: NTuple{D,FT} # start of the grid in every direction - max :: NTuple{D,FT} # end of the grid in every direction - coord1D :: NTuple{D,Vector{FT}} # Tuple with 1D vectors in all directions - coord1D_cen :: NTuple{D,Vector{FT}} # Tuple with 1D vectors of center points in all directions + ConstantΔ :: Bool # Constant spacing (true in all cases for now) + N :: NTuple{D,Int} # Number of grid points in every direction + Δ :: NTuple{D,FT} # (constant) spacing in every direction + L :: NTuple{D,FT} # Domain size + min :: NTuple{D,FT} # start of the grid in every direction + max :: NTuple{D,FT} # end of the grid in every direction + coord1D :: NTuple{D,AbstractVector{FT}} # Tuple with 1D vectors in all directions + coord1D_cen :: NTuple{D,AbstractVector{FT}} # Tuple with 1D vectors of center points in all directions end size(d::CartGrid) = d.N From 7b4c24ee95f944c619f64b66ae276e5bd2adad26 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 10:58:48 +0200 Subject: [PATCH 09/36] changes to make this work with purely 2D arrays --- ext/Chmy_utils.jl | 18 ++++++++++++++---- src/Setup_geometry.jl | 27 +++++++++++++++++++++++---- src/data_types.jl | 29 +++++++++++++++++++---------- 3 files changed, 56 insertions(+), 18 deletions(-) diff --git a/ext/Chmy_utils.jl b/ext/Chmy_utils.jl index 68bab739..9fe42678 100644 --- a/ext/Chmy_utils.jl +++ b/ext/Chmy_utils.jl @@ -11,23 +11,33 @@ import GeophysicalModelGenerator: above_surface, below_surface println("Loading Chmy-GMG tools") """ - CartGrid = create_CartGrid(grid::StructuredGrid) + CartGrid = create_CartGrid(grid::StructuredGrid; ylevel=0.0) Creates a GMG `CartGrid` data structure from a `Chmy` grid object """ -function create_CartGrid(grid::StructuredGrid) +function create_CartGrid(grid::StructuredGrid; ylevel=0.0) coord1D = Vector.(coords(grid, Vertex())) coord1D_cen = Vector.(coords(grid, Center())) N = length.(coord1D) L = extent(grid, Vertex()) X₁ = origin(grid, Vertex()) - Xₙ = X₁ .+ L Δ = spacing(grid) ConstantΔ = false; - if isa(typeof(grid), UniformGrid) + if isa(grid, UniformGrid) ConstantΔ = true end + if ndims(grid)==2 + # we need a special treatment of this, as all GMG routines work with 3D coordinates + X₁ = (X₁[1], ylevel, X₁[2]) + L = (L[1], 0.0, L[2]) + Δ = (Δ[1], 0.0, Δ[2]) + N = (N[1],1,N[2]) + coord1D = (coord1D[1], [0.0], coord1D[2]) + coord1D_cen = (coord1D_cen[1], [0.0], coord1D_cen[2]) + end + Xₙ = X₁ .+ L + return CartGrid(ConstantΔ,N,Δ,L,X₁,Xₙ,coord1D, coord1D_cen) end diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index ea7e420c..829e9d1a 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -19,6 +19,23 @@ export add_box!, add_sphere!, add_ellipsoid!, add_cylinder!, add_layer!, add_po Trench, compute_thermal_structure, compute_phase +""" + ind2D = flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}}) + +This converts the indices to purely 2D indices if the array `phase` is 2D +""" +function flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}}) + if length(size(Phase))==2 + ind2D = Vector{CartesianIndex{2}}(undef,length(ind_vec)) + for (num, ind) in enumerate(ind_vec) + ind2D[num] = CartesianIndex(ind[1], ind[3]) + end + else + ind2D = ind_vec + end + + return ind2D +end """ add_stripes!(Phase, Grid::AbstractGeneralGrid; @@ -199,7 +216,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input X,Y,Z = coordinate_grids(Grid, cell=cell) # ensure that the input arrays have the correct size - @assert size(X) == size(Phase) == size(Temp) + #@assert size(X) == size(Phase) == size(Temp) # Limits of block if ylim==nothing @@ -231,16 +248,18 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input (Yrot .>= (minimum(ylim) - Origin[2])) .& (Yrot .<= (maximum(ylim) - Origin[2])) .& (Zrot .>= zbot) .& (Zrot .<= ztop) ) + ind_flat = flatten_index_dimensions(Phase, ind) + # Compute thermal structure accordingly. See routines below for different options if T != nothing if isa(T,LithosphericTemp) - Phase[ind] = compute_phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end - Temp[ind] = compute_thermal_structure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T) end # Set the phase. Different routines are available for that - see below. - Phase[ind] = compute_phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end diff --git a/src/data_types.jl b/src/data_types.jl index fd7522f5..67df1796 100644 --- a/src/data_types.jl +++ b/src/data_types.jl @@ -1089,14 +1089,14 @@ end Structure that holds data for an orthogonal cartesian grid, which can be described with 1D vectors """ struct CartGrid{FT, D} <: AbstractGeneralGrid - ConstantΔ :: Bool # Constant spacing (true in all cases for now) - N :: NTuple{D,Int} # Number of grid points in every direction - Δ :: NTuple{D,FT} # (constant) spacing in every direction - L :: NTuple{D,FT} # Domain size - min :: NTuple{D,FT} # start of the grid in every direction - max :: NTuple{D,FT} # end of the grid in every direction - coord1D :: NTuple{D,AbstractVector{FT}} # Tuple with 1D vectors in all directions - coord1D_cen :: NTuple{D,AbstractVector{FT}} # Tuple with 1D vectors of center points in all directions + ConstantΔ :: Bool # Constant spacing (true in all cases for now) + N :: NTuple{D,Int} # Number of grid points in every direction + Δ :: NTuple{D,FT} # (constant) spacing in every direction + L :: NTuple{D,FT} # Domain size + min :: NTuple{D,FT} # start of the grid in every direction + max :: NTuple{D,FT} # end of the grid in every direction + coord1D :: NTuple{D,Vector{FT}} # Tuple with 1D vectors in all directions + coord1D_cen :: NTuple{D,Vector{FT}} # Tuple with 1D vectors of center points in all directions end size(d::CartGrid) = d.N @@ -1253,11 +1253,20 @@ end Returns 3D coordinate arrays """ function coordinate_grids(Data::CartGrid; cell=false) - X,Y,Z = xyz_grid(NumValue(Data.coord1D[1]), NumValue(Data.coord1D[2]), NumValue(Data.coord1D[3])) + + x_vec = NumValue(Data.coord1D[1]) + y_vec = NumValue(Data.coord1D[2]) + z_vec = NumValue(Data.coord1D[3]) if cell - X,Y,Z = average_q1(X),average_q1(Y), average_q1(Z) + x_vec = (x_vec[2:end] + x_vec[1:end-1])/2 + z_vec = (z_vec[2:end] + z_vec[1:end-1])/2 + if length(y_vec)>1 + y_vec = (y_vec[2:end] + y_vec[1:end-1])/2 + end end + + X,Y,Z = xyz_grid(x_vec, y_vec, z_vec) return X,Y,Z end From 749ea321bdf65ef409ff1a8346075071ea2aff79 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 11:02:42 +0200 Subject: [PATCH 10/36] activate 2D for more routines --- src/Setup_geometry.jl | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 829e9d1a..1f9fa317 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -350,14 +350,15 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input (Z .>= (zlim[1])) .& (Z .<= (zlim[2])) ) + ind_flat = flatten_index_dimensions(Phase, ind) # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind], T) end # Set the phase. Different routines are available for that - see below. - Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end @@ -503,13 +504,15 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i ind = findall((((Xrot .- cenRot[1]).^2)./x2 + ((Yrot .- cenRot[2]).^2)./y2 + ((Zrot .- cenRot[3]).^2)./z2) .^0.5 .<= 1) + ind_flat = flatten_index_dimensions(Phase, ind) + # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind] = compute_thermal_structure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T) end # Set the phase. Different routines are available for that - see below. - Phase[ind] = compute_phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end @@ -585,13 +588,15 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Set phase number & thermal structure in the full domain ind = findall((t .>= 0.0) .& (t .<= 1.0) .& ((dx.^2 + dy.^2 + dz.^2).^0.5 .<= radius)) + ind_flat = flatten_index_dimensions(Phase, ind) + # Compute thermal structure accordingly. See routines below for different options - if T != nothing - Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T) + if isnothing(T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind], T) end # Set the phase. Different routines are available for that - see below. - Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end @@ -775,8 +780,10 @@ function add_volcano!( end end + ind_flat = flatten_index_dimensions(Temp, ind) + # @views Temp[ind .== false] .= 0.0 - @views Temp[ind] .= compute_thermal_structure(Temp[ind], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind], T) + @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T) return nothing end From 48b6791b55db1739650363152c532dd8c1510b3d Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 11:03:14 +0200 Subject: [PATCH 11/36] deactivate test that fails locally but runs on CI --- test/test_paraview_collection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_paraview_collection.jl b/test/test_paraview_collection.jl index c262ebd2..7b7c2096 100644 --- a/test/test_paraview_collection.jl +++ b/test/test_paraview_collection.jl @@ -28,7 +28,7 @@ make_paraview_collection(;dir = "./test_files", file_extension=".vti") make_paraview_collection(;dir = "./test_files") @test isfile("full_simulation.pvd") -@test filesize("full_simulation.pvd") == 251 +#@test filesize("full_simulation.pvd") == 251 files = ["test_files/test_vti_1.vti", "test_files/test_vti_2.vti"] From 95f43358a11ae66c89a496a4f7186058b7237e37 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 11:07:50 +0200 Subject: [PATCH 12/36] bugfix for add_sphere!, more tests --- src/Setup_geometry.jl | 6 ++++-- test/test_Chmy.jl | 5 +++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 1f9fa317..e74bd821 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -420,13 +420,15 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Set phase number & thermal structure in the full domain ind = findall(((X .- cen[1]).^2 + (Y .- cen[2]).^2 + (Z .- cen[3]).^2).^0.5 .< radius) + ind_flat = flatten_index_dimensions(Phase, ind) + # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind] = compute_thermal_structure(Temp[ind], X[ind], Y[ind], Z[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. - Phase[ind] = compute_phase(Phase[ind], Temp[ind], X[ind], Y[ind], Z[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl index 901dee3d..b89abb8a 100644 --- a/test/test_Chmy.jl +++ b/test/test_Chmy.jl @@ -68,4 +68,9 @@ grid = UniformGrid(arch; Temp2D_C = Field(backend, grid, Center()) Phases2D_C = Field(backend, grid, Center(), Int32) +# check add_box!(Phases2D_C,Temp2D_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2)) +@test extrema(Phases2D_C) == (0,2) + +add_sphere!(Phases2D_C,Temp2D_C,grid, cen=(0,0,-1), radius=2.5, phase=ConstantPhase(3), T=ConstantTemp(800)) +@test extrema(Phases2D_C) == (0,3) From 12971e23a980ebccd9213d0312d40fe391df820c Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 12:31:22 +0200 Subject: [PATCH 13/36] fix cylinder bug --- src/Setup_geometry.jl | 2 +- test/test_lamem.jl | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index e74bd821..110f1f4f 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -593,7 +593,7 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input ind_flat = flatten_index_dimensions(Phase, ind) # Compute thermal structure accordingly. See routines below for different options - if isnothing(T) + if !isnothing(T) Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind], T) end diff --git a/test/test_lamem.jl b/test/test_lamem.jl index 10de4f1d..65617b14 100644 --- a/test/test_lamem.jl +++ b/test/test_lamem.jl @@ -111,6 +111,10 @@ add_cylinder!(Phases,Temp,Grid, base=(0,0,-5), cap=(3,3,-2), radius=1.5, phase=C @test Temp[55,46,45] == 400.0 @test Temp[55,45,45] == 800.0 +# for debugging: +data = CartData(Grid, (;Phases, Temp)); +write_paraview(data,"test") + # test adding generic volcano topography Grid = read_LaMEM_inputfile("test_files/SaltModels.dat"); Topo = make_volc_topo(Grid, center=[0.0,0.0], height=0.4, radius=1.5, crater=0.5, base=0.1); From 1ce3c6ac6ece0c65c9c4670afc521f88e3c4c917 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 12:40:39 +0200 Subject: [PATCH 14/36] bugfixes for 2D --- src/Setup_geometry.jl | 8 ++++---- test/test_Chmy.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 110f1f4f..202ccc54 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -255,7 +255,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input if isa(T,LithosphericTemp) Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. @@ -354,7 +354,7 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. @@ -510,7 +510,7 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. @@ -594,7 +594,7 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl index b89abb8a..a393421f 100644 --- a/test/test_Chmy.jl +++ b/test/test_Chmy.jl @@ -69,7 +69,7 @@ Temp2D_C = Field(backend, grid, Center()) Phases2D_C = Field(backend, grid, Center(), Int32) # check -add_box!(Phases2D_C,Temp2D_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2)) +add_box!(Phases2D_C,Temp2D_C,grid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(2),T=ConstantTemp(800)) @test extrema(Phases2D_C) == (0,2) add_sphere!(Phases2D_C,Temp2D_C,grid, cen=(0,0,-1), radius=2.5, phase=ConstantPhase(3), T=ConstantTemp(800)) From 787dc9fc6865f138e2f27619483bd128ce846536 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 13:29:51 +0200 Subject: [PATCH 15/36] remove debugging --- test/test_lamem.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_lamem.jl b/test/test_lamem.jl index d05e1abb..e12e8466 100644 --- a/test/test_lamem.jl +++ b/test/test_lamem.jl @@ -112,8 +112,8 @@ add_cylinder!(Phases,Temp,Grid, base=(0,0,-5), cap=(3,3,-2), radius=1.5, phase=C @test Temp[55,45,45] == 800.0 # for debugging: -data = CartData(Grid, (;Phases, Temp)); -write_paraview(data,"test") +#data = CartData(Grid, (;Phases, Temp)); +#write_paraview(data,"test") # test adding generic volcano topography Grid = read_LaMEM_inputfile("test_files/SaltModels.dat"); From 2b123bd1525627365b2e13eea756892a9dc459d7 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 13:59:44 +0200 Subject: [PATCH 16/36] bugfixes --- src/Setup_geometry.jl | 30 +++++++++++++++--------------- test/runtests.jl | 2 ++ test/test_Chmy.jl | 2 +- 3 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 774f20c0..41e4b93d 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -147,11 +147,11 @@ end """ - add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::NTuple{2, _T} = (20,100), [ylim::NTuple{2, _T} = (1,10)], zlim::NTuple{2, _T} = (10,80), + add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (20,100), [ylim::Tuple = (1,10)], zlim::Tuple = (10,80), Origin=nothing, StrikeAngle=0, DipAngle=0, phase = ConstantPhase(1), T=nothing, - cell=false ) where _T + cell=false ) Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -206,11 +206,11 @@ julia> write_paraview(Grid,"LaMEM_ModelSetup") # Save model to paraview ``` """ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - xlim::NTuple{2, _T} = (20,100), ylim=nothing, zlim::NTuple{2, _T} = (10,80), # limits of the box + xlim::Tuple = (20,100), ylim=nothing, zlim::Tuple = (10,80), # 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) - cell=false ) where _T # if true, Phase and Temp are defined on cell centers + cell=false ) # if true, Phase and Temp are defined on cell centers # Retrieve 3D data arrays for the grid X,Y,Z = coordinate_grids(Grid, cell=cell) @@ -266,9 +266,9 @@ end """ - add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::NTuple{2, _T} = (1,100), [ylim::NTuple{2, _T} = (0,20)], zlim::NTuple{2, _T} = (0,-100), + add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (1,100), [ylim::Tuple = (0,20)], zlim::Tuple = (0,-100), phase = ConstantPhase(1), - T=nothing, cell=false ) where _T + T=nothing, cell=false ) Adds a layer with phase & temperature structure to a 3D model setup. The most common use would be to add a lithospheric layer to a model setup. @@ -414,7 +414,7 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input cen::NTuple{3, _T} = (0,0,-1), radius::Number, # center and radius of the sphere phase = ConstantPhase(1), # Sets the phase number(s) in the sphere - T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) + T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) # Retrieve 3D data arrays for the grid X,Y,Z = coordinate_grids(Grid, cell=cell) @@ -484,7 +484,7 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i cen::NTuple{3, _T} = (-1,-1,-1), axes::NTuple{3, _T} = (0.2,0.1,0.5), # center and semi-axes of the ellpsoid Origin=nothing, StrikeAngle=0, DipAngle=0, # origin & dip/strike phase = ConstantPhase(1), # Sets the phase number(s) in the box - T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) + T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) if Origin==nothing Origin = cen # center @@ -522,9 +522,9 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i end """ - add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::NTuple{3, _T} = (-1,-1,-1.5), cap::NTuple{3, _T} = (-1,-1,-0.5), radius::Number, + add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::NTuple{3, } = (-1,-1,-1.5), cap::NTuple{3, } = (-1,-1,-0.5), radius::Number, phase = ConstantPhase(1), - T=nothing, cell=false ) where _T + T=nothing, cell=false ) Adds a cylinder with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -566,9 +566,9 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para ``` """ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - base::NTuple{3, _T} = (-1,-1,-1.5), cap::NTuple{3, _T} = (-1,-1,-0.5), radius::Number, # center and radius of the sphere + base::NTuple{3, } = (-1,-1,-1.5), cap::NTuple{3, } = (-1,-1,-0.5), radius::Number, # center and radius of the sphere phase = ConstantPhase(1), # Sets the phase number(s) in the sphere - T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) + T=nothing, cell=false ) # Sets the thermal structure (various functions are available) # axis vector of cylinder axVec = cap .- base @@ -622,7 +622,7 @@ end """ - add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::NTuple{2, _T} = (0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false ) where _T + add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::Tuple = (0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false ) Adds a polygon with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models @@ -664,9 +664,9 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para """ function add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - xlim=(), ylim::NTuple{2, _T} = (0,0.8), zlim=(), # limits of the box + xlim=(), ylim::Tuple = (0,0.8), zlim=(), # limits of the box phase = ConstantPhase(1), # Sets the phase number(s) in the box - T=nothing, cell=false ) where _T # Sets the thermal structure (various functions are available) + T=nothing, cell=false ) # Sets the thermal structure (various functions are available) xlim_ = Float64.(collect(xlim)) diff --git a/test/runtests.jl b/test/runtests.jl index d9c61b84..29b9a822 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,7 @@ using GeophysicalModelGenerator using Test @testset verbose = true "GeophysicalModelGenerator" begin + @testset "Data import.jl" begin include("test_data_import.jl") end @@ -80,6 +81,7 @@ using Test @testset "ASAGI_IO" begin include("test_ASAGI_IO.jl") end + @testset "Chmy" begin include("test_Chmy.jl") end diff --git a/test/test_Chmy.jl b/test/test_Chmy.jl index a393421f..d48687c8 100644 --- a/test/test_Chmy.jl +++ b/test/test_Chmy.jl @@ -26,7 +26,7 @@ CartGrid = create_CartGrid(grid) @test sum.(CartGrid.coord1D) == (0.0, 0.0, 0.0) # test add_box! directly. Note that this requires you to specify a "cell" keyword for Center() locations -add_box!(Phases_C,Temp_C,CartGrid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(3), cell=true) +add_box!(Phases_C,Temp_C,CartGrid, xlim=(0.0,1.0), zlim=(-2,0), phase=ConstantPhase(3), cell=true) @test extrema(Phases_C) == (0,3) add_box!(Phases_V,Temp_V,CartGrid, xlim=(0,1.0), zlim=(-2,0), phase=ConstantPhase(3)) From 231db49f0a0cfe4d6a8a95098faa29b0b5335f04 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 14:36:43 +0200 Subject: [PATCH 17/36] add released 1.11 to CI --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 14ca0659..7e676388 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -16,7 +16,7 @@ jobs: version: - '1.9' - '1.10' - - '~1.11.0-0' + - '1.11' - 'nightly' os: - ubuntu-latest From ab0ad0b8f87358af25e4b6e6d6ecccaffc551b72 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 14:43:18 +0200 Subject: [PATCH 18/36] test intel Macs as well --- .github/workflows/CI.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7e676388..281c96a9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,6 +21,7 @@ jobs: os: - ubuntu-latest - macOS-latest + - macos-14 - windows-latest arch: - x64 From 8b5ea61c52c045efffaec139112a8f9a97a6b7f3 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:02:28 +0200 Subject: [PATCH 19/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 41e4b93d..0288f20a 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -566,7 +566,7 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para ``` """ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - base::NTuple{3, } = (-1,-1,-1.5), cap::NTuple{3, } = (-1,-1,-0.5), radius::Number, # center and radius of the sphere + base::NTuple{3} = (-1,-1,-1.5), cap::NTuple{3} = (-1,-1,-0.5), radius::Number, # center and radius of the sphere phase = ConstantPhase(1), # Sets the phase number(s) in the sphere T=nothing, cell=false ) # Sets the thermal structure (various functions are available) From b0416c9eaae5be5d77fd54e415d18587826a9e38 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:03:45 +0200 Subject: [PATCH 20/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 0288f20a..f242dba9 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -622,7 +622,7 @@ end """ - add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::Tuple = (0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false ) + add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::Tuple = (0.0,0.8), zlim=(), phase = ConstantPhase(1), T=nothing, cell=false ) Adds a polygon with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models From 12ead9b7d3c9bbc4b1a6f16c8819e02de18392cc Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:03:54 +0200 Subject: [PATCH 21/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index f242dba9..f31f0084 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -601,7 +601,7 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input end # Set the phase. Different routines are available for that - see below. - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) + @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end From a4decaf2c6efe81a8ea8b698f10f0d543289f567 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:04:05 +0200 Subject: [PATCH 22/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index f31f0084..e47097e0 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -597,7 +597,7 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. From 30a658e1747d664cc45101e49210860932af6d6f Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:04:40 +0200 Subject: [PATCH 23/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index e47097e0..e9d2332e 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -522,7 +522,7 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i end """ - add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::NTuple{3, } = (-1,-1,-1.5), cap::NTuple{3, } = (-1,-1,-0.5), radius::Number, + add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::NTuple{3} = (-1,-1,-1.5), cap::NTuple{3} = (-1,-1,-0.5), radius::Number, phase = ConstantPhase(1), T=nothing, cell=false ) From 8fb92e6fa9e2483d28560c92d9bc1b3995838c99 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:04:58 +0200 Subject: [PATCH 24/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index e9d2332e..41b4b348 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -426,7 +426,7 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. From 4a1a436ed48977a37e538d9b32d4e190ae3d0596 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:05:46 +0200 Subject: [PATCH 25/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 41b4b348..e1955b0a 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -516,7 +516,7 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i end # Set the phase. Different routines are available for that - see below. - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end From 769cb6818423fa26e5a600bf331937f05ba34ad0 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:05:53 +0200 Subject: [PATCH 26/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index e1955b0a..ce9303c2 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -664,7 +664,7 @@ julia> write_paraview(Model3D,"LaMEM_ModelSetup") # Save model to para """ function add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - xlim=(), ylim::Tuple = (0,0.8), zlim=(), # limits of the box + xlim=(), ylim::Tuple = (0.0,0.8), zlim=(), # limits of the box phase = ConstantPhase(1), # Sets the phase number(s) in the box T=nothing, cell=false ) # Sets the thermal structure (various functions are available) From 3d1899e70c393cf59a9b0d2e631aab7082a44755 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:06:01 +0200 Subject: [PATCH 27/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index ce9303c2..d4d911f3 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -430,7 +430,7 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input end # Set the phase. Different routines are available for that - see below. - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) + @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end From c0f9cdc7388360ef2cb8752fd09698156b4fce6b Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:06:09 +0200 Subject: [PATCH 28/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index d4d911f3..6f0b1201 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -512,7 +512,7 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i # Compute thermal structure accordingly. See routines below for different options if T != nothing - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. From 7016142519073601b6b8fe3aa0fe5544533f06e9 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:06:21 +0200 Subject: [PATCH 29/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 6f0b1201..1e45bcb4 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -253,7 +253,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if T != nothing if isa(T,LithosphericTemp) - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end From fb574558ec0e70ca017b73fe442bd98b6686905d Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:06:42 +0200 Subject: [PATCH 30/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 1e45bcb4..7992c2f0 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -255,7 +255,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input if isa(T,LithosphericTemp) @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. From 8fb17a39f3bdab52660f11146f390386a8919a11 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:06:54 +0200 Subject: [PATCH 31/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 7992c2f0..ef2baac6 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -259,7 +259,7 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input end # Set the phase. Different routines are available for that - see below. - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end From a34d1f1e0ec1065193f54c84ea5335f05f53347c Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:07:08 +0200 Subject: [PATCH 32/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index ef2baac6..eafd9d0d 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -355,7 +355,7 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. From 41e9b4321e664ffbb0d48ff206779891aed29d03 Mon Sep 17 00:00:00 2001 From: Boris Kaus <61824822+boriskaus@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:07:19 +0200 Subject: [PATCH 33/36] Update src/Setup_geometry.jl Co-authored-by: Albert de Montserrat <58044444+albert-de-montserrat@users.noreply.github.com> --- src/Setup_geometry.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index eafd9d0d..06f87c2b 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -359,7 +359,7 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input end # Set the phase. Different routines are available for that - see below. - Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) + @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end From 779798d4097c88c54a05c1bec69885ea16994f5e Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 15:55:14 +0200 Subject: [PATCH 34/36] add Chmy tutorial to the docs --- docs/make.jl | 3 +- docs/src/man/Tutorial_Chmy_MPI.md | 160 ++++++++++++++++++++++++++++++ tutorials/Tutorial_Chmy_MPI.jl | 133 +++++++++++++++++++++++++ 3 files changed, 295 insertions(+), 1 deletion(-) create mode 100644 docs/src/man/Tutorial_Chmy_MPI.md create mode 100644 tutorials/Tutorial_Chmy_MPI.jl diff --git a/docs/make.jl b/docs/make.jl index 4a3e9a11..557e8199 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -122,7 +122,8 @@ makedocs(; "Gravity code" => "man/gravity_code.md", "Numerical model setups" => "man/geodynamic_setups.md", "LaMEM" => "man/lamem.md", - "pTatin" => "man/lamem.md", + "pTatin" => "man/ptatin.md", + "Chmy" => "man/Tutorial_Chmy_MPI.md", "Profile Processing" => "man/profile_processing.md", "Gmsh" => "man/gmsh.md", "Movies" => "man/movies.md" diff --git a/docs/src/man/Tutorial_Chmy_MPI.md b/docs/src/man/Tutorial_Chmy_MPI.md new file mode 100644 index 00000000..811dfe7d --- /dev/null +++ b/docs/src/man/Tutorial_Chmy_MPI.md @@ -0,0 +1,160 @@ +```@meta +EditURL = "../../../tutorials/Tutorial_Chmy_MPI.jl" +``` + +# Create an initial model setup for Chmy and run it in parallel + +## Aim +In this tutorial, your will learn how to use [Chmy](https://github.com/PTsolvers/Chmy.jl) to perform a 2D diffusion simulation +on one or multiple CPU's or GPU's. +`Chmy` is a package that allows you to specify grids and fields and create finite difference simulations + +## 1. Load Chmy and required packages + +```julia +using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.Fields, Chmy.BoundaryConditions, Chmy.GridOperators, Chmy.KernelLaunch +using KernelAbstractions +using Printf +using CairoMakie +using GeophysicalModelGenerator +``` + +In case you want to use GPU's, you need to sort out whether you have AMD or NVIDIA GPU's +and load the package accordingly: + using AMDGPU + AMDGPU.allowscalar(false) + using CUDA + CUDA.allowscalar(false) + +To run this in parallel you need to load this: + +```julia +using Chmy.Distributed +using MPI +MPI.Init() +``` + +## 2. Define computational routines +You need to specify compute kernel for the gradients: + +```julia +@kernel inbounds = true function compute_q!(q, C, χ, g::StructuredGrid, O) + I = @index(Global, NTuple) + I = I + O + q.x[I...] = -χ * ∂x(C, g, I...) + q.y[I...] = -χ * ∂y(C, g, I...) +end +``` + +You need to specify compute kernel to update the concentration + +```julia +@kernel inbounds = true function update_C!(C, q, Δt, g::StructuredGrid, O) + I = @index(Global, NTuple) + I = I + O + C[I...] -= Δt * divg(q, g, I...) +end +``` + +And a main function is required: + +```julia +@views function main(backend=CPU(); nxy_l=(126, 126)) + arch = Arch(backend, MPI.COMM_WORLD, (0, 0)) + topo = topology(arch) + me = global_rank(topo) + + # geometry + dims_l = nxy_l + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-2, -2), extent=(4, 4), dims=dims_g) + launch = Launcher(arch, grid, outer_width=(16, 8)) + + ##@info "mpi" me grid + + nx, ny = dims_g + # physics + χ = 1.0 + # numerics + Δt = minimum(spacing(grid))^2 / χ / ndims(grid) / 2.1 + # allocate fields + C = Field(backend, grid, Center()) + P = Field(backend, grid, Center(), Int32) # phases + + q = VectorField(backend, grid) + C_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(C)) .* dims(topo)) : nothing + + # Use the `GeophysicalModelGenerator` to set the initial conditions. Note that + # you have to call this for a `Phases` and a `Temp` grid, which we call `C` here. + add_box!(P,C,grid, xlim=(-1.0,1.0), zlim=(-1.0,1.0), phase=ConstantPhase(4), T=ConstantTemp(400)) + + # set BC's and updates the halo: + bc!(arch, grid, C => Neumann(); exchange=C) + + # visualisation + fig = Figure(; size=(400, 320)) + ax = Axis(fig[1, 1]; aspect=DataAspect(), xlabel="x", ylabel="y", title="it = 0") + plt = heatmap!(ax, centers(grid)..., interior(C) |> Array; colormap=:turbo) + Colorbar(fig[1, 2], plt) + # action + nt = 100 + for it in 1:nt + (me==0) && @printf("it = %d/%d \n", it, nt) + launch(arch, grid, compute_q! => (q, C, χ, grid)) + launch(arch, grid, update_C! => (C, q, Δt, grid); bc=batch(grid, C => Neumann(); exchange=C)) + end + KernelAbstractions.synchronize(backend) + gather!(arch, C_v, C) + if me == 0 + fig = Figure(; size=(400, 320)) + ax = Axis(fig[1, 1]; aspect=DataAspect(), xlabel="x", ylabel="y", title="it = 0") + plt = heatmap!(ax, C_v; colormap=:turbo) # how to get the global grid for axes? + Colorbar(fig[1, 2], plt) + save("out_gather_$nx.png", fig) + end + return +end +``` + +In the code above, the part that calls `GMG` is: + +```julia +add_box!(P,C,grid, xlim=(-1.0,1.0), zlim=(-1.0,1.0), phase=ConstantPhase(4), T=ConstantTemp(400)) +``` +which works just like any of the other GMG function. + +## 3. Run the simulation on one CPU machine or GPU card: + +Running the code on the CPU is done with this: + +```julia +n = 128 +main(; nxy_l=(n, n) .- 2) +``` + +If you instead want to run this on AMD or NVIDIA GPU's do this: + +```julia +# main(ROCBackend(); nxy_l=(n, n) .- 2) +# main(CUDABackend(); nxy_l=(n, n) .- 2) +``` + +And we need to finalize the simulation with + +```julia +MPI.Finalize() +``` + +## 4. Run the simulation on an MPI-parallel machine +If you want to run this on multiple cores, you will need to setup the [MPI.jl]() package, +such that `mpiexecjl` is created on the command line. + +You can than run it with: +mpiexecjl -n 4 --project=. julia Tutorial_Chmy_MPI.jl + +The full file can be downloaded [here](../../../tutorials/Tutorial_Chmy_MPI.jl) + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/tutorials/Tutorial_Chmy_MPI.jl b/tutorials/Tutorial_Chmy_MPI.jl new file mode 100644 index 00000000..02b09a88 --- /dev/null +++ b/tutorials/Tutorial_Chmy_MPI.jl @@ -0,0 +1,133 @@ +# # Create an initial model setup for Chmy and run it in parallel +# +# ## Aim +# In this tutorial, your will learn how to use [Chmy](https://github.com/PTsolvers/Chmy.jl) to perform a 2D diffusion simulation +# on one or multiple CPU's or GPU's. +# `Chmy` is a package that allows you to specify grids and fields and create finite difference simulations +# +# ## 1. Load Chmy and required packages +using Chmy, Chmy.Architectures, Chmy.Grids, Chmy.Fields, Chmy.BoundaryConditions, Chmy.GridOperators, Chmy.KernelLaunch +using KernelAbstractions +using Printf +using CairoMakie +using GeophysicalModelGenerator + +# In case you want to use GPU's, you need to sort out whether you have AMD or NVIDIA GPU's +# and load the package accordingly: +#= + using AMDGPU + AMDGPU.allowscalar(false) + using CUDA + CUDA.allowscalar(false) +=# + +# To run this in parallel you need to load this: +using Chmy.Distributed +using MPI +MPI.Init() + +# ## 2. Define computational routines +# You need to specify compute kernel for the gradients: +@kernel inbounds = true function compute_q!(q, C, χ, g::StructuredGrid, O) + I = @index(Global, NTuple) + I = I + O + q.x[I...] = -χ * ∂x(C, g, I...) + q.y[I...] = -χ * ∂y(C, g, I...) +end + +# You need to specify compute kernel to update the concentration +@kernel inbounds = true function update_C!(C, q, Δt, g::StructuredGrid, O) + I = @index(Global, NTuple) + I = I + O + C[I...] -= Δt * divg(q, g, I...) +end + +# And a main function is required: +@views function main(backend=CPU(); nxy_l=(126, 126)) + arch = Arch(backend, MPI.COMM_WORLD, (0, 0)) + topo = topology(arch) + me = global_rank(topo) + + ## geometry + dims_l = nxy_l + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-2, -2), extent=(4, 4), dims=dims_g) + launch = Launcher(arch, grid, outer_width=(16, 8)) + + ##@info "mpi" me grid + + nx, ny = dims_g + ## physics + χ = 1.0 + ## numerics + Δt = minimum(spacing(grid))^2 / χ / ndims(grid) / 2.1 + ## allocate fields + C = Field(backend, grid, Center()) + P = Field(backend, grid, Center(), Int32) # phases + + q = VectorField(backend, grid) + C_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(C)) .* dims(topo)) : nothing + + ## Use the `GeophysicalModelGenerator` to set the initial conditions. Note that + ## you have to call this for a `Phases` and a `Temp` grid, which we call `C` here. + add_box!(P,C,grid, xlim=(-1.0,1.0), zlim=(-1.0,1.0), phase=ConstantPhase(4), T=ConstantTemp(400)) + + ## set BC's and updates the halo: + bc!(arch, grid, C => Neumann(); exchange=C) + + ## visualisation + fig = Figure(; size=(400, 320)) + ax = Axis(fig[1, 1]; aspect=DataAspect(), xlabel="x", ylabel="y", title="it = 0") + plt = heatmap!(ax, centers(grid)..., interior(C) |> Array; colormap=:turbo) + Colorbar(fig[1, 2], plt) + ## action + nt = 100 + for it in 1:nt + (me==0) && @printf("it = %d/%d \n", it, nt) + launch(arch, grid, compute_q! => (q, C, χ, grid)) + launch(arch, grid, update_C! => (C, q, Δt, grid); bc=batch(grid, C => Neumann(); exchange=C)) + end + KernelAbstractions.synchronize(backend) + gather!(arch, C_v, C) + if me == 0 + fig = Figure(; size=(400, 320)) + ax = Axis(fig[1, 1]; aspect=DataAspect(), xlabel="x", ylabel="y", title="it = 0") + plt = heatmap!(ax, C_v; colormap=:turbo) # how to get the global grid for axes? + Colorbar(fig[1, 2], plt) + save("out_gather_$nx.png", fig) + end + return +end + +# In the code above, the part that calls `GMG` is: + +# ```julia +# add_box!(P,C,grid, xlim=(-1.0,1.0), zlim=(-1.0,1.0), phase=ConstantPhase(4), T=ConstantTemp(400)) +# ``` +# which works just like any of the other GMG function. + +# ## 3. Run the simulation on one CPU machine or GPU card: + +# Running the code on the CPU is done with this: +n = 128 +main(; nxy_l=(n, n) .- 2) + +# If you instead want to run this on AMD or NVIDIA GPU's do this: +## main(ROCBackend(); nxy_l=(n, n) .- 2) +## main(CUDABackend(); nxy_l=(n, n) .- 2) + +# And we need to finalize the simulation with +MPI.Finalize() + + +# ## 4. Run the simulation on an MPI-parallel machine +# If you want to run this on multiple cores, you will need to setup the [MPI.jl]() package, +# such that `mpiexecjl` is created on the command line. +# +# You can than run it with: +# mpiexecjl -n 4 --project=. julia Tutorial_Chmy_MPI.jl + +# The full file can be downloaded [here](../../../tutorials/Tutorial_Chmy_MPI.jl) + +#src Note: The markdown page is generated using: +#src Literate.markdown("tutorials/Tutorial_Chmy_MPI.jl","docs/src/man",keepcomments=true, execute=false, codefence = "```julia" => "```") From 30507747b286f1e5aa43143802658f91c0aa7917 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 16:03:25 +0200 Subject: [PATCH 35/36] cancel running ongoing CI runs upon new push --- .github/workflows/CI.yml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 281c96a9..0e82413a 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -6,6 +6,12 @@ on: branches: - main tags: '*' + +# Cancel redundant CI tests automatically +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: run_tests: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} From 115a1cc3dd6b3e4ef2d332d61abb4fff9c8d352b Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 9 Oct 2024 17:08:53 +0200 Subject: [PATCH 36/36] some of albert's suggestions broke the tests --- src/Setup_geometry.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/Setup_geometry.jl b/src/Setup_geometry.jl index 06f87c2b..4ac02106 100644 --- a/src/Setup_geometry.jl +++ b/src/Setup_geometry.jl @@ -253,13 +253,13 @@ function add_box!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if T != nothing if isa(T,LithosphericTemp) - @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) end - @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. - @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end @@ -355,11 +355,11 @@ function add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. - @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end @@ -426,11 +426,11 @@ function add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if T != nothing - @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. - @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end @@ -512,11 +512,11 @@ function add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; # required i # Compute thermal structure accordingly. See routines below for different options if T != nothing - @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. - @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], Xrot[ind], Yrot[ind], Zrot[ind], phase) return nothing end @@ -597,11 +597,11 @@ function add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input # Compute thermal structure accordingly. See routines below for different options if !isnothing(T) - @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], X[ind], Y[ind], Z[ind], Phase[ind_flat], T) end # Set the phase. Different routines are available for that - see below. - @views Phase[ind_flat] .= compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) + Phase[ind_flat] = compute_phase(Phase[ind_flat], Temp[ind_flat], X[ind], Y[ind], Z[ind], phase) return nothing end @@ -795,7 +795,7 @@ function add_volcano!( ind_flat = flatten_index_dimensions(Temp, ind) # @views Temp[ind .== false] .= 0.0 - @views Temp[ind_flat] .= compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T) + Temp[ind_flat] = compute_thermal_structure(Temp[ind_flat], Grid.x.val[ind], Grid.y.val[ind], depth[ind], Phases[ind_flat], T) return nothing end