Skip to content

Commit

Permalink
update scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
@aelligp
aelligp committed Sep 24, 2024
1 parent eb3de4c commit 1ba0067
Show file tree
Hide file tree
Showing 5 changed files with 300 additions and 428 deletions.
122 changes: 13 additions & 109 deletions SillModelSetup.jl
View file
Original file line number Diff line number Diff line change
@@ -1,63 +1,35 @@
## Model Setup
using GeophysicalModelGenerator

function SillSetup(nx, ny, nz; sticky_air=5)
Lx = Ly = 50
function SillSetup(nx, ny, nz)
Lx = Ly = 300
x = range(0.0, Lx, nx)
y = range(0.0, Ly, 2)
z = range(-25, sticky_air, nz)
z = range(-250, 0.0, nz)
Grid = CartData(xyz_grid(x, y, z))

# 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(4, nx, 2, nz)
Phases = fill(1, nx, 2, nz)

# In many (geodynamic) models, one also has to define the temperature, so lets define it as well
Temp = fill(0.0, nx, 2, nz)
Temp = fill(600+273, nx, 2, nz)

add_box!(
Phases,
Temp,
Grid;
xlim=(minimum(Grid.x.val), maximum(Grid.x.val)),
ylim=(minimum(Grid.y.val), maximum(Grid.y.val)),
zlim=(minimum(Grid.z.val), 0.0),
phase=LithosphericPhases(; Layers=[30], Phases=[1]),
T=HalfspaceCoolingTemp(; Age=20),
)

# add_volcano!(Phases, Temp, Grid;
# volcanic_phase = 1,
# center = (mean(Grid.x.val), 0.0),
# height = 3,
# radius = 5,
# crater = 0.5,
# base = 0.0,
# background = nothing,
# T = HalfspaceCoolingTemp(Age=20)
# )

add_ellipsoid!(
Phases,
Temp,
Grid;
cen=(mean(Grid.x.val), 0, -5.0),
axes=(3.5, 2.5, 2.0),
zlim=(-175, -75.0),
phase=ConstantPhase(2),
T=ConstantTemp(; T=1000),
T = ConstantTemp(; T=1000+273),
)

# add_sphere!(Phases, Temp, Grid;
# cen = (mean(Grid.x.val), 0,-5.0),
# radius = 0.5,
# cen = (mean(Grid.x.val), 0,-90.0),
# radius = 5,
# phase = ConstantPhase(3),
# T = ConstantTemp(T=1100)
# )
# add_cylinder!(Phases, Temp, Grid;
# base = (mean(Grid.x.val), 0, -3.25),
# cap = (mean(Grid.x.val), 0, 3.00),
# radius = 0.05,
# phase = ConstantPhase(2),
# # T = LinearTemp(Ttop=20, Tbot=1000),
# # T = ConstantTemp(T=800),
# T = ConstantTemp(T=1000),
# T = ConstantTemp(T=1100+273)
# )

Grid = addfield(Grid, (; Phases, Temp))
Expand All @@ -67,74 +39,6 @@ function SillSetup(nx, ny, nz; sticky_air=5)

ph = Phases[:, 1, :]
T = Temp[:, 1, :]
# write_paraview(Grid, "Volcano2D")
return li, origin, ph, T, Grid
end

function simple_setup_no_FS2D(nx, ny, nz)
Lx = Ly = 40

x = range(-Lx / 2, Lx / 2, nx)
z = range(-20, 0, nz)
Grid = CartData(xyz_grid(x, 0, z))

# 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(4, nx, 1, nz)

# In many (geodynamic) models, one also has to define the temperature, so lets define it as well
Temp = fill(20.0, nx, 1, nz)

add_box!(
Phases,
Temp,
Grid;
xlim=(minimum(Grid.x.val), maximum(Grid.x.val)),
zlim=(minimum(Grid.z.val), maximum(Grid.z.val)),
phase=LithosphericPhases(; Layers=[30], Phases=[1]),
T=HalfspaceCoolingTemp(; Age=20),
)

add_ellipsoid!(
Phases,
Temp,
Grid;
cen=(0, 0, -5),
axes=(2.5, 2.5, 2.5 / 2),
phase=ConstantPhase(2),
T=ConstantTemp(; T=1000),
)
add_sphere!(
Phases,
Temp,
Grid;
cen=(0, 0, -5.0),
radius=0.5,
phase=ConstantPhase(3),
T=ConstantTemp(; T=1200),
)
add_cylinder!(
Phases,
Temp,
Grid;
base=(0, 0, -4.0),
cap=(0, 0, 4.0),
radius=0.5,
phase=ConstantPhase(2),
T=LinearTemp(; Ttop=20, Tbot=1000),
)

surf = Grid.z.val .> 0.0
Temp[surf] .= 20.0
Phases[surf] .= 4
@. Temp[Phases == 4] = 20
@. Temp = max(Temp, 20)
Grid = addfield(Grid, (; Phases, Temp))

li = (abs(last(x) - first(x)), abs(last(z) - first(z))) .* 1e3
origin = (x[1], z[1]) .* 1e3

ph = Phases
T = Temp
write_paraview(Grid, "Volcano2D")
return li, origin, ph, T
return li, origin, ph, T, Grid
end
108 changes: 108 additions & 0 deletions SillRheology.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@
# from "Fingerprinting secondary mantle plumes", Cloetingh et al. 2022

function init_rheologies(CharDim)
linear_viscosity_rhy = LinearMeltViscosity(A = -8.1590, B = 2.4050e+04K, T0 = -430.9606K)#,η0=1e0Pa*s)
linear_viscosity_bas = LinearMeltViscosity(A = -9.6012, B = 1.3374e+04K, T0 = 307.8043K)#, η0=1e3Pa*s)

# Define rheolgy struct
rheology = (
# Name = "Rhyolite",
SetMaterialParams(;
Phase = 1,
Density = MeltDependent_Density(ρsolid=ConstantDensity(ρ=2700kg / m^3),ρmelt=T_Density(ρ0=2300kg / m^3)),
HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity(), Q_L=400e3J/kg),
Conductivity = ConstantConductivity(),
CompositeRheology = CompositeRheology((linear_viscosity_rhy,)),
Melting = MeltingParam_Smooth3rdOrder(a=3043.0,b=-10552.0,c=12204.9,d=-4709.0),
CharDim = CharDim,
),
# Name = "Basaltic_Sill",
SetMaterialParams(;
Phase = 2,
Density = MeltDependent_Density(ρsolid=ConstantDensity(ρ=3000kg/m^3),ρmelt=T_Density(ρ0=2800kg / m^3)),
HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity(), Q_L=400e3J/kg),
Conductivity = ConstantConductivity(),
CompositeRheology = CompositeRheology((linear_viscosity_bas,)),
Melting = MeltingParam_Smooth3rdOrder(),
CharDim = CharDim,
),
# # Name = "Basaltic_Sill_thermal_Anomaly",
# SetMaterialParams(;
# Phase = 3,
# Density = MeltDependent_Density(ρsolid=ConstantDensity(ρ=2900kg/m^3),ρmelt=T_Density(ρ0=2800kg / m^3)),
# HeatCapacity = Latent_HeatCapacity(Cp=ConstantHeatCapacity(), Q_L=400e3J/kg),
# Conductivity = ConstantConductivity(),
# CompositeRheology = CompositeRheology((linear_viscosity_bas,)),
# Melting = MeltingParam_Smooth3rdOrder(),
# CharDim = CharDim,
# ),
)
end


# function init_phases!(phases, particles)
# ni = size(phases)

# @parallel_indices (i, j) function init_phases!(phases, px, py, index)
# @inbounds for ip in JustPIC._2D.cellaxes(phases)
# # quick escape
# @index(index[ip, i, j]) == 0 && continue

# x = @index px[ip, i, j]
# depth = -(@index py[ip, i, j])
# @index phases[ip, i, j] = 1.0

# if 0.1e3 < depth ≤ 0.2e3
# @index phases[ip, i, j] = 2.0

# end
# end
# return nothing
# end

# @parallel (@idx ni) init_phases!(phases, particles.coords..., particles.index)
# end
function init_phases!(phases, phase_grid, particles, xvi)
ni = size(phases)
@parallel (@idx ni) _init_phases!(
phases, phase_grid, particles.coords, particles.index, xvi
)
end

@parallel_indices (I...) function _init_phases!(
phases, phase_grid, pcoords::NTuple{N,T}, index, xvi
) where {N,T}
ni = size(phases)

for ip in cellaxes(phases)
# quick escape
@index(index[ip, I...]) == 0 && continue

pᵢ = ntuple(Val(N)) do i
@index pcoords[i][ip, I...]
end

d = Inf # distance to the nearest particle
particle_phase = -1
for offi in 0:1, offj in 0:1
ii = I[1] + offi
jj = I[2] + offj

!(ii ≤ ni[1]) && continue
!(jj ≤ ni[2]) && continue

xvᵢ = (xvi[1][ii], xvi[2][jj])
d_ijk = √(sum((pᵢ[i] - xvᵢ[i])^2 for i in 1:N))
if d_ijk < d
d = d_ijk
particle_phase = phase_grid[ii, jj]
end
# if pᵢ[end] > 0.0 && phase_grid[ii, jj] > 1.0
# particle_phase = 4.0
# end
end
@index phases[ip, I...] = Float64(particle_phase)
end

return nothing
end
Loading

0 comments on commit 1ba0067

Please sign in to comment.