From 57b5d099f9aedbd748c0b5ab746dd3d40271ebdd Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Wed, 30 Oct 2024 09:55:43 +0100 Subject: [PATCH] progress push --- .../diffusion/diffusion2D_multiphase_MPI.jl | 10 -- miniapps/subduction/3D/Subduction3D_MPI.jl | 170 +++++++++++------- miniapps/subduction/3D/Subduction3D_setup.jl | 18 +- 3 files changed, 115 insertions(+), 83 deletions(-) diff --git a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase_MPI.jl b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase_MPI.jl index c03e5a36..947ba253 100644 --- a/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase_MPI.jl +++ b/miniapps/benchmarks/thermal_diffusion/diffusion/diffusion2D_multiphase_MPI.jl @@ -73,16 +73,6 @@ function diffusion_2D(figdir; nx=32, ny=32, lx=100e3, ly=100e3, Cp0=1.2e3, K0=3. ttot = 1 * Myr # total simulation time dt = 50 * kyr # physical time step - # init_mpi = JustRelax.MPI.Initialized() ? false : true - # igg = IGG(init_global_grid(nx, ny, 1; select_device=false, init_MPI = init_mpi)...) - - # # Physical domain - # ni = (nx, ny) - # li = (lx, ly) # domain length in x- and y- - # di = @. li / (nx_g(), ny_g()) # grid step in x- and -y - # grid = Geometry(ni, li; origin = (0, -ly)) - # (; xci, xvi) = grid # nodes at the center and vertices of the cells - # Physical domain ni = nx, ny li = lx, ly # domain length in x- and y- diff --git a/miniapps/subduction/3D/Subduction3D_MPI.jl b/miniapps/subduction/3D/Subduction3D_MPI.jl index b0388d69..b889f4b1 100644 --- a/miniapps/subduction/3D/Subduction3D_MPI.jl +++ b/miniapps/subduction/3D/Subduction3D_MPI.jl @@ -39,7 +39,7 @@ function main3D(li, origin, phases_GMG, igg; nx=16, ny=16, nz=16, figdir="figs3D # Physical domain ------------------------------------ ni = nx, ny, nz # number of cells - di = @. li / ni # grid steps + di = @. li / (nx_g(), ny_g(), nz_g()) # grid steps grid = Geometry(ni, li; origin = origin) (; xci, xvi) = grid # nodes at the center and vertices of the cells # ---------------------------------------------------- @@ -175,9 +175,11 @@ function main3D(li, origin, phases_GMG, igg; nx=16, ny=16, nz=16, figdir="figs3D ) ); end - println("Stokes solver time ") - println(" Total time: $t_stokes s") - println(" Time/iteration: $(t_stokes / out.iter) s") + if igg.me == 0 + println("Stokes solver time ") + println(" Total time: $t_stokes s") + println(" Time/iteration: $(t_stokes / out.iter) s") + end tensor_invariant!(stokes.ε) dt = compute_dt(stokes, di) # ------------------------------ @@ -219,76 +221,101 @@ function main3D(li, origin, phases_GMG, igg; nx=16, ny=16, nz=16, figdir="figs3D inject_particles_phase!(particles, pPhases, (), (), xvi) # update phase ratios update_phase_ratios!(phase_ratios, particles, xci, xvi, pPhases) + if igg.me == 0 + @show it += 1 + t += dt + end + if do_vtk - @show it += 1 - t += dt - - # Data I/O and plotting --------------------- - if igg.me == 0 && (it == 1 || rem(it, 20) == 0) - # checkpointing(figdir, stokes, thermal.T, η, t) - - if do_vtk - velocity2vertex!(Vx_v, Vy_v, Vz_v, @velocity(stokes)...) - phase_vertex = [argmax(p) for p in Array(phase_ratios.vertex)] - phase_center = [argmax(p) for p in Array(phase_ratios.center)] - - # MPI - #centers - @views P_nohalo .= Array(stokes.P[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views τII_nohalo .= Array(stokes.τ.II[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views η_vep_nohalo .= Array(stokes.viscosity.η_vep[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views εII_nohalo .= Array(stokes.ε.II[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views phases_c_v .= Array(phase_center[2:end-1, 2:end-1, 2:end-1]) - gather!(P_nohalo, P_v) - gather!(τII_nohalo, τII_v) - gather!(η_vep_nohalo, η_vep_v) - gather!(εII_nohalo, εII_v) - gather!(phases_c_nohalo, phases_c_v) - #vertices - @views Vxv_nohalo .= Array(Vx_v[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views Vyv_nohalo .= Array(Vy_v[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views Vzv_nohalo .= Array(Vz_v[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views T_nohalo .= Array(thermal.T[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo - @views phase_v_nohalo .= Array(phase_vertex[2:end-1, 2:end-1, 2:end-1]) - gather!(Vxv_nohalo, Vxv_v) - gather!(Vyv_nohalo, Vyv_v) - gather!(Vzv_nohalo, Vzv_v) - gather!(T_nohalo, T_v) - gather!(phase_v_nohalo, phases_v_v) - - - - - - + velocity2vertex!(Vx_v, Vy_v, Vz_v, @velocity(stokes)...) data_v = (; - T = T_v, - phases_v = phases_v_v, + T = zeros(ni.+1...), + phase_vertex = [argmax(p) for p in Array(phase_ratios.center)] ) data_c = (; - P = P_v, - τII = τII_v, - εII = εII_v, - η = η_vep_v, - phases = phases_c_v - - + # P = Array(stokes.P), + # τII = Array(stokes.τ.II), + # εII = Array(stokes.ε.II), + η = Array(stokes.viscosity.η), + phase_center = [argmax(p) for p in Array(phase_ratios.center)] ) velocity_v = ( - Array(Vxv_v), - Array(Vyv_v), - Array(Vzv_v), + Array(Vx_v), + Array(Vy_v), + Array(Vz_v), ) save_vtk( - joinpath(vtk_dir, "vtk_" * lpad("$it", 6, "0")), - xvi_v, - xci_v, + joinpath(vtk_dir, "vtk_" * lpad("$(it)_$(igg.me)", 6, "0")), + xvi, + xci, data_v, data_c, velocity_v ) - end end + # Data I/O and plotting --------------------- + # if igg.me == 0 && (it == 1 || rem(it, 20) == 0) + # # checkpointing(figdir, stokes, thermal.T, η, t) + + # if do_vtk + # velocity2vertex!(Vx_v, Vy_v, Vz_v, @velocity(stokes)...) + # phase_vertex = [argmax(p) for p in Array(phase_ratios.vertex)] + # phase_center = [argmax(p) for p in Array(phase_ratios.center)] + # println("pv= $(size(phase_vertex)), pc= $(size(phase_center))") + # println("stokes P = $(size(stokes.P)), T = $(size(thermal.T))") + # # MPI + # #centers + # @views P_nohalo .= Array(stokes.P[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # @views τII_nohalo .= Array(stokes.τ.II[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # @views η_vep_nohalo .= Array(stokes.viscosity.η_vep[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # @views εII_nohalo .= Array(stokes.ε.II[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # # @views phases_c_v .= Array(phase_center[2:end-1, 2:end-1, 2:end-1]) + # gather!(P_nohalo, P_v) + # gather!(τII_nohalo, τII_v) + # gather!(η_vep_nohalo, η_vep_v) + # gather!(εII_nohalo, εII_v) + # # gather!(phases_c_nohalo, phases_c_v) + # #vertices + # @views Vxv_nohalo .= Array(Vx_v[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # @views Vyv_nohalo .= Array(Vy_v[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # @views Vzv_nohalo .= Array(Vz_v[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # @views T_nohalo .= Array(thermal.T[2:end-1, 2:end-1, 2:end-1]) # Copy data to CPU removing the halo + # # @views phase_v_nohalo .= Array(phase_vertex[2:end-1, 2:end-1, 2:end-1]) + # gather!(Vxv_nohalo, Vxv_v) + # gather!(Vyv_nohalo, Vyv_v) + # gather!(Vzv_nohalo, Vzv_v) + # gather!(T_nohalo, T_v) + # # gather!(phase_v_nohalo, phases_v_v) + + + # # data_v = (; + # # T = T_v, + # # # phases_v = phases_v_v, + # # ) + # # data_c = (; + # # P = P_v, + # # τII = τII_v, + # # εII = εII_v, + # # η = η_vep_v, + # # # phases = phases_c_v + + + # # ) + # # velocity_v = ( + # # Array(Vxv_v), + # # Array(Vyv_v), + # # Array(Vzv_v), + # # ) + # # save_vtk( + # # joinpath(vtk_dir, "vtk_" * lpad("$it", 6, "0")), + # # xvi_v, + # # xci_v, + # # data_v, + # # data_c, + # # velocity_v + # # ) + # end + # end # ------------------------------ end @@ -301,13 +328,28 @@ do_vtk = true # set to true to generate VTK files for ParaView # nx,ny,nz = 50, 50, 50 nx,ny,nz = 32,32,32 # nx,ny,nz = 128, 32, 64 -li, origin, phases_GMG, = GMG_only(nx+1, ny+1, nz+1) igg = if !(JustRelax.MPI.Initialized()) # initialize (or not) MPI grid IGG(init_global_grid(nx, ny, nz; init_MPI= true)...) else igg end +# Physical domain ------------------------------------ +x = range(-3960, 500, nx_g()); +y = range(0, 2640, ny_g()); +air_thickness = 0.0 +z = range(-660, air_thickness, nz_g()); +origin = (x[1], y[1], z[1]) +li = (abs(last(x)-first(x)), abs(last(y)-first(y)), abs(last(z)-first(z))) + +ni = nx, ny, nz # number of cells +di = @. li / (nx_g(), ny_g(), nz_g()) # grid steps +grid = Geometry(ni, li; origin = origin) +(; xci, xvi) = grid # nodes at the center and vertices of the cells +li, origin, phases_GMG, = GMG_only(xvi,nx, ny, nz) +# li, origin, phases_GMG, = GMG_only(nx_g()+1, ny_g()+1, nz_g()+1) +# li, origin, phases_GMG, = GMG_only(Grid, nx+1, ny+1, nz+1) + # (Path)/folder where output data and figures are stored -figdir = "Subduction3D_$(nx)x$(ny)x$(nz)" -main3D(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, nz = nz, do_vtk = do_vtk); +# figdir = "Subduction3D_$(nx_g())x$(ny_g())x$(nz_g())" +# main3D(li, origin, phases_GMG, igg; figdir = figdir, nx = nx, ny = ny, nz = nz, do_vtk = do_vtk); diff --git a/miniapps/subduction/3D/Subduction3D_setup.jl b/miniapps/subduction/3D/Subduction3D_setup.jl index 28323bca..b4954e38 100644 --- a/miniapps/subduction/3D/Subduction3D_setup.jl +++ b/miniapps/subduction/3D/Subduction3D_setup.jl @@ -1,17 +1,17 @@ using GeophysicalModelGenerator -function GMG_only(nx, ny, nz) +function GMG_only(xvi, nx, ny, nz) # nx,ny,nz = 99, 33, 66 # nx, ny, nz = (nx,ny,nz) .+ 1 - x = range(-3960, 500, nx); - y = range(0, 2640, ny); - air_thickness = 0.0 - z = range(-660, air_thickness, nz); + x = range(minimum(xvi[1]), maximum(xvi[1]), nx); + y = range(minimum(xvi[2]), maximum(xvi[2]), ny); + z = range(minimum(xvi[3]), maximum(xvi[3]), nz); Grid = CartData(xyz_grid(x,y,z)); + println("Grid_$(igg.me) created with $(nx) x $(ny) x $(nz) cells, $Grid") # Now we create an integer array that will hold the `Phases` information (which usually refers to the material or rock type in the simulation) Phases = fill(3, nx, ny, nz); - + # In many (geodynamic) models, one also has to define the temperature, so lets define it as well Temp = fill(1350.0, nx, ny, nz); @@ -28,7 +28,7 @@ function GMG_only(nx, ny, nz) # # And an an inclined part: # lith = LithosphericPhases(Layers=[200], Phases=[1 2], Tlab=1250) - # add_box!(Phases, Temp, Grid; xlim=(0,300).-1000, ylim=(0, 1000.0), zlim=(-80.0, 10.0), phase = lith, + # add_box!(Phases, Temp, Grid; xlim=(0,300).-1000, ylim=(0, 1000.0), zlim=(-80.0, 10.0), phase = lith, # # Origin=(-1000,0,0), # T=SpreadingRateTemp(SpreadingVel=0, MORside="left"), DipAngle=20, StrikeAngle=0); @@ -39,9 +39,9 @@ function GMG_only(nx, ny, nz) Grid = addfield(Grid,(;Phases, Temp)) # Which looks like - write_paraview(Grid,"Initial_Setup_Subduction"); + write_paraview(Grid,"Initial_Setup_Subduction_$(igg.me)"); - surf = Grid.z.val .> 0.0 + surf = Grid.z.val .> 0.0 Phases[surf] .= 4 li = (abs(last(x)-first(x)), abs(last(y)-first(y)), abs(last(z)-first(z))) .* 1e3