Merge pull request #141 from JuliaGeodynamics/bk-chmy

Integration with Chmy.jl

Author: boriskaus

.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 }}
@@ -16,11 +22,12 @@ jobs:
         version:
           - '1.9'
           - '1.10'
-          - '~1.11.0-0'
+          - '1.11'
           - 'nightly'
         os:
           - ubuntu-latest
           - macOS-latest
+          - macos-14
           - windows-latest
         arch:
           - x64

Project.toml

@@ -1,7 +1,7 @@
 name = "GeophysicalModelGenerator"
 uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742"
 authors = ["Boris Kaus", "Marcel Thielmann"]
-version = "0.7.7"
+version = "0.7.8"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -32,23 +32,26 @@ 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"
 
 [compat]
+Chmy = "0.1.20"
 Colors = "0.9 - 0.12"
 Downloads = "1"
 FFMPEG = "0.4"
 FileIO = "1"
-GDAL_jll = "300.900.0 - 301.900.0"
-GLMakie = "0.10"
-GMT = "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"
@@ -63,15 +66,17 @@ NearestNeighbors = "0.2 - 0.4"
 Parameters = "0.9 - 0.12"
 SpecialFunctions = "1.0, 2"
 StaticArrays = "1"
-WhereTheWaterFlows = "0.10"
+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"]
+test = ["Test", "GMT", "StableRNGs", "GridapGmsh", "Chmy","KernelAbstractions"]

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"            

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).*
+

ext/Chmy_utils.jl

@@ -0,0 +1,108 @@
+#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, 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")
+
+"""
+    CartGrid = create_CartGrid(grid::StructuredGrid; ylevel=0.0)
+
+Creates a GMG `CartGrid` data structure from a `Chmy` grid object
+"""
+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())  
+    Δ           = spacing(grid)
+    ConstantΔ   = false;
+    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
+
+# 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
+
+
+# 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` which can be either on vertices 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
+
+
+end # end of module