From c0b25917e6c8d92d03be9606d925e6dc612c5c97 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Wed, 12 Jul 2023 12:46:47 +0200 Subject: [PATCH] - adds a function to create the initial model setup (also in parallel) - simplified GMG routines - uses `files` as default - add plotting routines that are loaded when GLMakie is loaded --- src/LaMEM.jl | 14 +++- src/LaMEM_ModelGeneration/ErrorChecking.jl | 6 ++ src/LaMEM_ModelGeneration/GMG_interface.jl | 63 ++++-------------- .../LaMEM_ModelGeneration.jl | 7 +- src/LaMEM_ModelGeneration/Model.jl | 45 +++++++++++-- src/LaMEM_ModelGeneration/ModelSetup.jl | 2 +- src/LaMEM_ModelGeneration/Utils.jl | 34 +++++++++- src/MakieExt.jl | 65 +++++++++++++++++++ 8 files changed, 177 insertions(+), 59 deletions(-) create mode 100644 src/MakieExt.jl diff --git a/src/LaMEM.jl b/src/LaMEM.jl index c20aae5..67cd549 100644 --- a/src/LaMEM.jl +++ b/src/LaMEM.jl @@ -2,6 +2,7 @@ module LaMEM using GeoParams using .GeoParams +using Requires export NO_units, GEO_units, SI_units, km, m, Pa, Pas, kg, cm, yr #export GeoParams @@ -18,19 +19,28 @@ using .Run export run_lamem, run_lamem_save_grid, mpiexec export remove_popup_messages_mac, show_paths_LaMEM +# Functions that help running LaMEM directly from julia include("LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl") using .LaMEM_Model -export Model, Write_LaMEM_InputFile, +export Model, Write_LaMEM_InputFile, create_initialsetup, Scaling, Grid, Time, FreeSurface, BoundaryConditions, SolutionParams, Solver, ModelSetup, geom_Sphere, Output, Materials, Phase, Softening, PhaseTransition, Dike, - add_phase!, rm_phase!, rm_last_phase!, add_petsc!, add_softening!, + add_phase!, rm_phase!, rm_last_phase!, replace_phase!, add_petsc!, add_softening!, add_phasetransition!, add_dike!, add_geom! + using .Run.LaMEM_jll export LaMEM_jll # export LaMEM_jll as well & directories +function __init__() + @require GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin + @eval include("MakieExt.jl") + end +end + + end # module diff --git a/src/LaMEM_ModelGeneration/ErrorChecking.jl b/src/LaMEM_ModelGeneration/ErrorChecking.jl index 7a003e4..77c38c8 100644 --- a/src/LaMEM_ModelGeneration/ErrorChecking.jl +++ b/src/LaMEM_ModelGeneration/ErrorChecking.jl @@ -17,6 +17,12 @@ function Check_LaMEM_Model(m::Model) error("If you use internal geometries to set phases, you need to at least specify one internal geometry object. Example: add_geom!(model, geom_Sphere())") end + + if (m.ModelSetup.msetup=="files") && diff([extrema(m.Grid.Phases)...])[1]==0 && diff([extrema(m.Grid.Temp)...])[1]==0 + error("Your initial `Temp` grid is constant, as is your initial `Phases` grid. + You do need to set some variability to see action, for example with the GMG function `AddSphere!(model,cen=(0.0,0.0,0.0), radius=(0.15, ))` ") + end + return nothing end diff --git a/src/LaMEM_ModelGeneration/GMG_interface.jl b/src/LaMEM_ModelGeneration/GMG_interface.jl index e8aedc2..0606623 100644 --- a/src/LaMEM_ModelGeneration/GMG_interface.jl +++ b/src/LaMEM_ModelGeneration/GMG_interface.jl @@ -8,19 +8,13 @@ import GeophysicalModelGenerator: AddBox!, AddSphere!, AddEllipsoid!, AddCylinde AddBox!(model::Model; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2}, Origin=nothing, StrikeAngle=0, DipAngle=0, phase = ConstantPhase(1), - T=nothing ) + T=nothing ) Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models See the documentation of the GMG routine """ -function AddBox!(model::Model; xlim=Tuple{2}, ylim=Tuple{2}, zlim=Tuple{2}, - Origin=nothing, StrikeAngle=0, DipAngle=0, phase = ConstantPhase(1), T=nothing ) - -return AddBox!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; xlim=xlim, ylim=ylim, zlim=zlim, - Origin=Origin, StrikeAngle=StrikeAngle, DipAngle=DipAngle, phase = phase, T=T ) - -end +AddBox!(model::Model; kwargs...) = AddBox!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) """ @@ -29,58 +23,29 @@ end See the documentation of the GMG routine """ -function AddSphere!(model::Model; cen=Tuple{3}, radius=Tuple{1}, phase = ConstantPhase(1), T=nothing) - - AddSphere!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid, cen=cen, radius=radius, phase=phase, T=T) - - return nothing -end - - +AddSphere!(model::Model; kwargs...) = AddSphere!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) """ - AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1}, # center and radius of the sphere - phase = ConstantPhase(1), # Sets the phase number(s) in the sphere - T=nothing ) # Sets the thermal structure (various fucntions are available) + AddCylinder!(model::Model; # required input + base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1}, # center and radius of the sphere + phase = ConstantPhase(1), # Sets the phase number(s) in the sphere + T=nothing ) # Sets the thermal structure (various fucntions are available) See the documentation of the GMG routine """ -function AddCylinder!(Phase, Temp, Grid::AbstractGeneralGrid; # required input - base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1}, # center and radius of the sphere - phase = ConstantPhase(1), # Sets the phase number(s) in the sphere - T=nothing ) # Sets the thermal structure (various fucntions are available) - - - AddCylinder!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid, - base=base, cap=cap, radius=radius, # center and radius of the sphere - phase = phase, # Sets the phase number(s) in the sphere - T=T ) - - return nothing -end +AddCylinder!(model::Model; kwargs...) = AddCylinder!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) + """ - AddEllipsoid!(model::Model; # required input - cen=Tuple{3}, axes=Tuple{3}, # 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 + AddEllipsoid!(model::Model; # required input + cen=Tuple{3}, axes=Tuple{3}, # 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 ) See the documentation of the GMG routine """ -function AddEllipsoid!(model::Model; # required input - cen=Tuple{3}, axes=Tuple{3}, # 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 ) # Sets the thermal structure (various fucntions are available) - - AddEllipsoid!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid, - cen=cen, axes=axes, Origin=Origin, StrikeAngle=StrikeAngle, DipAngle=DipAngle, - phase=phase, T=T) - - return nothing -end \ No newline at end of file +AddEllipsoid!(model::Model; kwargs...) = AddEllipsoid!(model.Grid.Phases, model.Grid.Temp, model.Grid.Grid; kwargs...) diff --git a/src/LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl b/src/LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl index 5ed019e..023b5df 100644 --- a/src/LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl +++ b/src/LaMEM_ModelGeneration/LaMEM_ModelGeneration.jl @@ -83,10 +83,15 @@ include("Utils.jl") export add_phase!, rm_phase!, rm_last_phase!, add_petsc!, add_softening!, add_phasetransition!, add_dike!, add_geom! - include("ErrorChecking.jl") export Check_LaMEM_Model +# optional dependency +#if !isdefined(Base, :get_extension) +# include("MakieExt.jl") +#end + + end diff --git a/src/LaMEM_ModelGeneration/Model.jl b/src/LaMEM_ModelGeneration/Model.jl index 5833a4e..6081f5b 100644 --- a/src/LaMEM_ModelGeneration/Model.jl +++ b/src/LaMEM_ModelGeneration/Model.jl @@ -1,10 +1,11 @@ # This is the main LaMEM Model struct using GeophysicalModelGenerator.GeoParams -import LaMEM.Run: run_lamem +import LaMEM.Run: run_lamem, run_lamem_save_grid +using LaMEM.Run.LaMEM_jll -export Model, Write_LaMEM_InputFiles +export Model, Write_LaMEM_InputFiles, create_initialsetup -""" +"""; Model Structure that holds all the information to create a LaMEM input file @@ -176,15 +177,49 @@ Performs a LaMEM run for the parameters that are specified in `model` """ function run_lamem(model::Model, cores::Int64=1, args::String=""; wait=true) + create_initialsetup(model, cores, args); + + run_lamem(model.Output.param_file_name, cores, args; wait=wait) + + return nothing +end + +""" + create_initialsetup(model::Model, cores::Int64=1, args::String="") + +Creates the initial model setup of LaMEM from `model`, which includes: +- Writing the LaMEM (*.dat) input file +- Write the VTK file (if requested when `model.Output.write_VTK_setup=true`) +- Write the marker files to disk (if `model.ModelSetup.msetup="files"`) + +""" +function create_initialsetup(model::Model, cores::Int64=1, args::String="") + + # Move to the working directory + cur_dir = pwd() + if !isempty(model.Output.out_dir) + if !exist(model.Output.out_dir); mkdir(model.Output.out_dir); end # create directory if needed + cd(model.Output.out_dir) + end + Write_LaMEM_InputFile(model, model.Output.param_file_name) if model.ModelSetup.msetup=="files" # write marker files to disk before running LaMEM + Model3D = CartData(model.Grid.Grid, (Phases=model.Grid.Phases,Temp=model.Grid.Temp)); - end + if cores>1 + # PartFile = CreatePartitioningFile(model.Output.param_file_name, cores, LaMEM_options=args) + + PartFile = run_lamem_save_grid(model.Output.param_file_name, cores) - run_lamem(model.Output.param_file_name, cores, args; wait=wait) + Save_LaMEMMarkersParallel(Model3D, PartitioningFile=PartFile) + else + Save_LaMEMMarkersParallel(Model3D) + end + end + cd(cur_dir) return nothing end diff --git a/src/LaMEM_ModelGeneration/ModelSetup.jl b/src/LaMEM_ModelGeneration/ModelSetup.jl index b13b635..f9ed4bf 100644 --- a/src/LaMEM_ModelGeneration/ModelSetup.jl +++ b/src/LaMEM_ModelGeneration/ModelSetup.jl @@ -9,7 +9,7 @@ export ModelSetup, Write_LaMEM_InputFile, geom_Sphere, set_geom! """ Base.@kwdef mutable struct ModelSetup "Setup type - can be `geom` (phases are assigned from geometric primitives), `files` (from julia input), `polygons` (from geomIO input, which requires `poly_file` to be specified) " - msetup::String = "geom" + msetup::String = "files" "add random noise to the particle location" rand_noise::Int64 = 1 diff --git a/src/LaMEM_ModelGeneration/Utils.jl b/src/LaMEM_ModelGeneration/Utils.jl index 60d8eec..5ba8ab6 100644 --- a/src/LaMEM_ModelGeneration/Utils.jl +++ b/src/LaMEM_ModelGeneration/Utils.jl @@ -1,5 +1,6 @@ # Contains a number of useful functions -export add_phase!, rm_phase!, rm_last_phase!, add_softening!, add_phasetransition!, +export add_phase!, rm_phase!, rm_last_phase!, replace_phase!, + add_softening!, add_phasetransition!, add_dike!, add_geom! @@ -12,6 +13,17 @@ function add_phase!(model::Model, phase::Phase) return nothing end +""" + add_phase!(model::Model, phases...) +Add several phases @ once. +""" +function add_phase!(model::Model, phases...) + for phase in phases + push!(model.Materials.Phases, phase); + end +end + + """ rm_last_phase!(model::Model, phase::Phase) This removes the last added `phase` from `model` @@ -35,6 +47,26 @@ function rm_phase!(model::Model, ID::Int64) end +""" + replace_phase!(model::Model, phase_new::Phase; ID::Int64=nothing, Name::String=nothing) + +This replaces a `phase` within a LaMEM Model Setup `model` with `phase_new` either based on its `Name` or `ID`. +Note that it is expected that only one such phase is present in the current setup. +""" +function replace_phase!(model::Model, phase_new::Phase; ID::Union{Nothing,Int64}=nothing, Name::Union{Nothing,String}=nothing) + id_vec = [phase.ID for phase in model.Materials.Phases] + name_vec = [phase.Name for phase in model.Materials.Phases] + if !isnothing(ID) + id = findfirst(id_vec .== ID) + elseif !isnothing(Name) + id = findfirst(name_vec .== Name) + end + model.Materials.Phases[id] = phase_new + + return nothing +end + + """ add_petsc!(model::Model, option::String) diff --git a/src/MakieExt.jl b/src/MakieExt.jl new file mode 100644 index 0000000..c3e58e7 --- /dev/null +++ b/src/MakieExt.jl @@ -0,0 +1,65 @@ +# Plotting extensions, that are only loaded when GLMakie is available +println("adding GLMakie plotting extensions for LaMEM") + +using LaMEM, GeophysicalModelGenerator +using .GLMakie + +export plot_initial_setup + +""" + plot_initial_setup(model::LaMEM.Model; x=nothing, y=nothing, z=nothing, phases=true) + +Plots a 2D setup through the LaMEM model setup `model`. if `phases=true` we plot phases; otherwise temperature +""" +function plot_initial_setup(model::LaMEM.Model; x=nothing, y=nothing, z=nothing, phases=true) + + Model3D = CartData(model.Grid.Grid, (Phases=model.Grid.Phases,Temp=model.Grid.Temp)); + + if !isnothing(z); Depth_level = z*km; else Depth_level=nothing; end + if !isnothing(x); Lon_level = x; else Lon_level=nothing; end + if !isnothing(y); Lat_level = y; else Lat_level=nothing; end + + Cross = CrossSectionVolume(Model3D, Depth_level=Depth_level, Lat_level=Lat_level, Lon_level=Lon_level) + + # Flatten arrays + coordinates + if !isnothing(x) + Phase2D, Temp2D = dropdims(Cross.fields.Phases,dims=1), dropdims(Cross.fields.Temp,dims=1) + X2D, Z2D = dropdims(Cross.y.val,dims=1), dropdims(Cross.z.val,dims=1) + x_vec, z_vec = X2D[:,1], Z2D[1,:] + title_str = "x = $x" + x_str,z_str = "y","z" + elseif !isnothing(y) + Phase2D, Temp2D = dropdims(Cross.fields.Phases,dims=2), dropdims(Cross.fields.Temp,dims=2) + X2D, Z2D = dropdims(Cross.x.val,dims=2), dropdims(Cross.z.val,dims=2) + x_vec, z_vec = X2D[:,1], Z2D[1,:] + title_str = "y = $y" + x_str,z_str = "x","z" + elseif !isnothing(z) + Phase2D, Temp2D = dropdims(Cross.fields.Phases,dims=3), dropdims(Cross.fields.Temp,dims=3) + X2D, Z2D = dropdims(Cross.x.val,dims=3), dropdims(Cross.y.val,dims=3) + + x_vec, z_vec = X2D[:,1], Z2D[1,:] + title_str = "z = $z" + x_str,z_str = "x","y" + end + + + # Create figure + if phases + data = Phase2D + title_str = "Phases "*title_str + else + data = Temp2D + title_str = "Temperature "*title_str + end + + fig = Figure() + ax = Axis(fig[1, 1], aspect = DataAspect(), title = title_str, xlabel=x_str, ylabel=z_str) + hm = heatmap!(ax, x_vec, z_vec, data) + + Colorbar(fig[:, end+1], hm) + + display(fig) + + return nothing +end