Skip to content

Commit

Permalink
up
Browse files Browse the repository at this point in the history
albert-de-montserrat committed Jun 24, 2024
1 parent 3e949ad commit 791c108
Showing 3 changed files with 29 additions and 26 deletions.
14 changes: 8 additions & 6 deletions Volcano3D.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# const isCUDA = false
const isCUDA = true
const isCUDA = false
# const isCUDA = true

@static if isCUDA
using CUDA
@@ -126,7 +126,7 @@ function main3D(
)
εbg = nondimensionalize(1e-15 / s, CharDim)
pureshear_bc!(
stokes, xci, xvi, εbg, backend
stokes, xci, xvi, εbg, backend_JR
)

flow_bcs = FlowBoundaryConditions(;
@@ -316,7 +316,9 @@ function main3D(
εyy = Array(ustrip.(dimensionalize(stokes.ε.yy, s^-1,CharDim))),
εzz = Array(ustrip.(dimensionalize(stokes.ε.zz, s^-1,CharDim))),
εII = Array(ustrip.(dimensionalize(stokes.ε.II, s^-1,CharDim))),
η = Array(ustrip.(dimensionalize(stokes.viscosity.η_vep,Pa*s,CharDim))),
ρ = Array(ustrip.(dimensionalize(ρg[end] ./ rheology[1].Gravity[1].g.val,kg/m^3,CharDim))),
η = Array(ustrip.(dimensionalize(stokes.viscosity.η,Pa*s,CharDim))),
η_vep = Array(ustrip.(dimensionalize(stokes.viscosity.η_vep,Pa*s,CharDim))),
)
velocity_v = (
Array(ustrip.(dimensionalize(Vx_v,cm/yr,CharDim))),
@@ -361,7 +363,7 @@ end

figdir = "Volcano3D"
do_vtk = true # set to true to generate VTK files for ParaView
n = 50
n = 256
nx = n
ny = n
nz = n
@@ -373,4 +375,4 @@ else
end

# run main script
main3D(igg, li_GMG, origin_GMG, phases_GMG, T_GMG; figdir=figdir, nx=nx, ny=ny, nz=nz, do_vtk = do_vtk);
# main3D(igg, li_GMG, origin_GMG, phases_GMG, T_GMG; figdir=figdir, nx=nx, ny=ny, nz=nz, do_vtk = do_vtk);
2 changes: 1 addition & 1 deletion VolcanoRheology.jl
Original file line number Diff line number Diff line change
@@ -79,7 +79,7 @@ function init_rheology(CharDim; is_compressible = false)
Density = ConstantDensity= 100kg/m^3,),
# Density = ConstantDensity(ρ=1kg/m^3,),
HeatCapacity = ConstantHeatCapacity(; Cp=1000J / kg / K),
Conductivity = ConstantConductivity(; k=15Watt / K / m),
Conductivity = ConstantConductivity(; k=1e-1Watt / K / m),
LatentHeat = ConstantLatentHeat(; Q_L=0.0J / kg),
ShearHeat = ConstantShearheating(0.0NoUnits),
CompositeRheology = CompositeRheology((creep_air,)),
39 changes: 20 additions & 19 deletions VolcanoSetup.jl
Original file line number Diff line number Diff line change
@@ -3,31 +3,29 @@ using GeophysicalModelGenerator, GMT
using Interpolations

function Etna_topo(N)
Topo = import_topo([14.5, 15.5, 37.2, 38.2], file="@earth_relief_01m")
# import topo from GMT's server
Topo = import_topo([14.5, 15.5, 37.2, 38.2], file="@earth_relief_03s")
# extract 2D array with the topographic elevation
surf = Topo.depth.val[:,:,1]
# l = maximum(abs.(Topo.depth.val))

l = maximum(abs.(Topo.depth.val))
# heatmap(Topo.depth.val[:,:,1], colormap=:oleron, colorrange= (-l,l))
# GLMakie.surface(surf, colormap=:oleron, colorrange= (-l,l))

##
# our model is in [km] -> need to convert from lat-long
x = LinRange(14.5, 15.5, N) |> collect
y = LinRange(37.2, 38.2, N) |> collect

X = (
Topo.lon.val[:, 1, 1],
Topo.lat.val[1, :, 1],
)

# interpolate from GMT's surface to the resolution nx × ny resolution of our grid
itp = interpolate(X, surf, Gridded(Linear()))
surf_interp = [itp(x, y) for x in x, y in y]

# GLMakie.surface(surf_interp, colormap=:oleron, colorrange= (-l,l))

# compute the x and y size of our cartesian model
lat_dist = extrema(X[1]) |> collect |> diff |> first |> abs
long_dist = extrema(X[2]) |> collect |> diff |> first |> abs
Ly = lat_dist * 110.574
Lx = long_dist * 111.320*cos(lat_dist)
Lx = lat_dist * 110.574
Ly = long_dist * 111.320*cos(lat_dist)

return surf_interp, Lx, Ly
end
@@ -38,14 +36,14 @@ function volcano_setup(N)
nx = ny = nz = N
x = range(-Lx / 2, Lx / 2, nx);
y = range(-Ly / 2, Ly / 2, ny);
z = range(-50, 10, nz);
z = range(-40, 10, 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, ny, nz);
Phases = fill(1, 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);
Temp = fill(20.0, nx, ny, nz);

lith = LithosphericPhases(Layers=[15 45 100], Phases=[1 2 3])

@@ -56,11 +54,6 @@ function volcano_setup(N)
phase = lith,
T = HalfspaceCoolingTemp(Age=20)
)

id = falses(nx, ny, nz)
for k in axes(topo_volcano,3), j in axes(topo_volcano, 2), i in axes(topo_volcano,1)
id[i, j, k] = Grid.z.val[i, j, k] > topo_volcano[i, j]
end

# add_volcano!(Phases, Temp, Grid;
# volcanic_phase = 1,
@@ -79,6 +72,14 @@ function volcano_setup(N)
phase = ConstantPhase(3),
)

# id = falses(nx, ny, nz)
for k in axes(Grid.z.val, 3), j in axes(topo_volcano, 2), i in axes(topo_volcano,1)
# id[i, j, k] = Grid.z.val[i, j, k] > topo_volcano[i, j]
if Grid.z.val[i, j, k] > topo_volcano[i, j]
Phases[i, j, k] = 4
end
end

@. Temp[Phases == 3] += 500
@. Temp[Phases == 4] = 20
@. Temp = max(Temp, 20)

0 comments on commit 791c108

Please sign in to comment.