Skip to content

Commit c84a751

Browse files
danieldouglas92danieldouglas92
committed
T-Model implementation
1 parent f35b5f7 commit c84a751

File tree

4 files changed

+355
-22
lines changed

4 files changed

+355
-22
lines changed

‎cookbooks/landlab/2D-T-Models/T-Model-LandLab.py‎

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,185 @@
1+
2+
from mpi4py import MPI
3+
import json
4+
import os
5+
import numpy as np
6+
import landlab
7+
from landlab.components import LinearDiffuser
8+
from landlab.components import FlowAccumulator
9+
from landlab.components import StreamPowerEroder
10+
11+
from landlab.io.legacy_vtk import write_legacy_vtk
12+
13+
current_time = 0
14+
15+
comm = None
16+
17+
model_grid = None
18+
elevation = None
19+
linear_diffuser = None
20+
flow_accumulator = None
21+
stream_power_eroder = None
22+
23+
s2yr = 365.25 * 24 * 3600
24+
timestep = 0
25+
vtks = []
26+
27+
def initialize(comm_handle):
28+
if not comm_handle is None:
29+
# Convert the handle back to an MPI communicator
30+
global comm
31+
comm = MPI.Comm.f2py(comm_handle)
32+
33+
rank = comm.Get_rank()
34+
size = comm.Get_size()
35+
36+
print(f"Python: Hello from Rank {rank} of {size}")
37+
38+
data = 1
39+
globalsum = comm.allreduce(data, op=MPI.SUM)
40+
if comm.rank == 0:
41+
print(f"\tPython: testing communication; sum {globalsum}")
42+
else:
43+
print("Python: running sequentially!")
44+
45+
def finalize():
46+
pass
47+
48+
49+
# Run the Landlab simulation from the current time to end_time and return
50+
# the new topographic elevation (in m) at each local node.
51+
# dict_variable_name_to_value_in_nodes is a dictionary mapping variables
52+
# (x velocity, y velocity, temperature, etc.) to an array of values in each
53+
# node.
54+
def update_until(end_time, dict_variable_name_to_value_in_nodes):
55+
global current_time, elevation, linear_diffuser, flow_accumulator, stream_power_eroder, timestep
56+
dt = end_time - current_time
57+
timestep += 1
58+
59+
deposition_erosion = np.zeros(model_grid.number_of_nodes)
60+
61+
x_velocity = dict_variable_name_to_value_in_nodes["x velocity"]
62+
y_velocity = dict_variable_name_to_value_in_nodes["y velocity"]
63+
z_velocity = dict_variable_name_to_value_in_nodes["z velocity"]
64+
65+
# In the current setup, the y_velocity vector is ONLY nonzero along y=0. We need
66+
# to project the velocity outwards along all fixed values of y, assigning the same
67+
# velocity to all nodes that share the same x value.
68+
69+
# This code below first extract the velocities along y=0 (where they are correctly coming
70+
# from ASPECT), then it finds all of the unique y values in the LandLab mesh, and then assigned
71+
# the y=0 velocities to these y-values. This will only work for structured RasterGrids.
72+
73+
slice_y_velocity = y_velocity[model_grid.y_of_node == 0]
74+
for y in np.unique(model_grid.y_of_node):
75+
y_velocity[model_grid.y_of_node == y] = slice_y_velocity
76+
vertical_velocity = y_velocity
77+
78+
# Substepping for surface processes
79+
if dt>0:
80+
n_substeps = 10
81+
sub_dt = dt / n_substeps
82+
for _ in range(n_substeps):
83+
84+
# TODO:
85+
elevation_before = elevation.copy()
86+
87+
flow_accumulator.run_one_step()
88+
stream_power_eroder.run_one_step(sub_dt)
89+
linear_diffuser.run_one_step(sub_dt)
90+
91+
elevation += vertical_velocity * sub_dt
92+
93+
deposition_erosion += elevation - elevation_before
94+
pass
95+
96+
filename = f"./output-2D-T-MODELS/landlab_{str(timestep).zfill(3)}.vtk"
97+
print("Writing output VTK file...", filename)
98+
vtk_file = write_legacy_vtk(path=filename, grid=model_grid, clobber=True, z_at_node=elevation)
99+
100+
vtks.append((current_time, filename))
101+
102+
if True:
103+
# write vtk.series file (ParaView supports legacy VTK in this format)
104+
with open("./output-2D-T-MODELS/landlab.vtk.series", "w") as f:
105+
series = {
106+
"file-series-version": "1.0",
107+
"files": [
108+
{"name": os.path.basename(filename), "time": time / s2yr}
109+
for time, filename in vtks
110+
]
111+
}
112+
json.dump(series, f, indent=2)
113+
114+
current_time = end_time
115+
print("Max elevation:", np.max(elevation), "Min elevation:", np.min(elevation))
116+
117+
# This setup does not do any averaging of the topography across the LandLab mesh before returninn
118+
# the change in elevation used to deform the ASPECT mesh. Averaging will be straight forward to
119+
# do going forward, but for the purposes of this test I want to make sure the ASPECT mesh is
120+
# deforming to the exact same topography as the LandLab mesh where they intersect.
121+
return deposition_erosion
122+
123+
def set_mesh_information(dict_grid_information):
124+
global model_grid, elevation
125+
126+
if not model_grid:
127+
print("* Creating RasterModelGrid ...")
128+
x_extent = dict_grid_information["Mesh X extent"]
129+
y_extent = dict_grid_information["Mesh Y extent"]
130+
spacing = dict_grid_information["Mesh Spacing"]
131+
132+
nrows = int(y_extent / spacing) + 1 # number of node rows
133+
ncols = int(x_extent / spacing) + 1 # number of node columns
134+
135+
# For 2D models, shift the LandLab mesh so that it is centered on the y-axis.
136+
model_grid = landlab.RasterModelGrid((nrows, ncols), xy_spacing=(spacing, spacing), xy_of_lower_left=(0, -y_extent / 2))
137+
138+
model_grid.set_closed_boundaries_at_grid_edges(right_is_closed=True,
139+
left_is_closed=False,
140+
top_is_closed=True,
141+
bottom_is_closed=True)
142+
143+
print("* Creating topographic elevation ...")
144+
# Initialize topography array with zeros
145+
elevation = model_grid.add_zeros("topographic__elevation", at="node")
146+
# Assign 1000 m high topography with random noise
147+
np.random.seed(0) # For reproducibility
148+
elevation += 1000 + np.random.rand(model_grid.number_of_nodes)
149+
150+
print("* Done")
151+
152+
# Return the x coordinates of the locally owned nodes on this
153+
# MPI rank. grid_id is always 0.
154+
def get_grid_x(grid_id):
155+
global model_grid
156+
return model_grid.node_x
157+
158+
# Return the y coordinates of the locally owned nodes on this
159+
# MPI rank. grid_id is always 0.
160+
def get_grid_y(grid_id):
161+
global model_grid
162+
return model_grid.node_y
163+
164+
def initialize_landlab_components(landlab_component_parameters):
165+
global model_grid, elevation, linear_diffuser, flow_accumulator, stream_power_eroder
166+
167+
D = landlab_component_parameters["Hillslope coefficient"]
168+
K_sp = landlab_component_parameters["Stream power erodibility coefficient"]
169+
flow_director = "D8" # Only flow director supported for HexGrids
170+
171+
172+
print("* Creating LinearDiffuser ... with D =", D)
173+
linear_diffuser = LinearDiffuser(model_grid, linear_diffusivity=D)
174+
print("* Creating FlowAccumulator ... with flow director =", flow_director)
175+
flow_accumulator = FlowAccumulator(model_grid, elevation, flow_director=flow_director)
176+
print("* Creating StreamPowerEroder ... with K_sp =", K_sp)
177+
stream_power_eroder = StreamPowerEroder(model_grid, K_sp=K_sp)
178+
179+
pass
180+
# Return the initial topography at the start of the simulation
181+
# in each node.
182+
def get_initial_topography(grid_id):
183+
global elevation
184+
return elevation
185+

‎cookbooks/landlab/2D-T-Models/simple_2D.prm‎

Lines changed: 138 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,138 @@
1+
# This parameter file sets up a model to test the interfacing between
2+
# the surface processes code LandLab and ASPECT. The ASPECT model is a 2D
3+
# vertical slice model, and the LandLab model is a 2D horizontal surface model.
4+
# This model tests that the interfacing between LandLab and ASPECT performs
5+
# correctly in this 2D, "T-Model" configuration. The ASPECT mesh and the LandLab
6+
# mesh are the same resoltuion, and since there is no averaging of topography
7+
# on the LandLab side, the topography should exactly match where the two models
8+
# intersect.
9+
set Dimension = 2
10+
11+
12+
set Use years instead of seconds = true
13+
set End time = 1000e3
14+
set Maximum time step = 10e3
15+
16+
17+
set Nonlinear solver scheme = single Advection, single Stokes
18+
set Pressure normalization = surface
19+
set Surface pressure = 0
20+
set Output directory = output-2D-T-MODELS
21+
22+
# 2.5 km mesh for ASPECT
23+
subsection Geometry model
24+
set Model name = box
25+
26+
subsection Box
27+
28+
set X extent = 100e3
29+
set Y extent = 100e3
30+
31+
set X repetitions = 40
32+
set Y repetitions = 40
33+
34+
end
35+
end
36+
37+
# Temperature effects are ignored
38+
subsection Initial temperature model
39+
set Model name = function
40+
subsection Function
41+
set Function expression = 273
42+
end
43+
end
44+
45+
subsection Boundary temperature model
46+
set Fixed temperature boundary indicators = bottom, top
47+
set List of model names = initial temperature
48+
end
49+
50+
# Free slip on all boundaries
51+
subsection Boundary velocity model
52+
set Prescribed velocity boundary indicators = left: function, right: function, bottom: function
53+
subsection Function
54+
set Coordinate system = cartesian
55+
set Variable names = x, y
56+
set Function constants = uplift=0.01
57+
set Function expression = 0; x * uplift / 100e3
58+
end
59+
end
60+
61+
subsection Mesh deformation
62+
set Mesh deformation boundary indicators = top: Landlab
63+
set Additional tangential mesh velocity boundary indicators = left, right
64+
65+
subsection Landlab
66+
set MPI ranks for Landlab = 1
67+
set Script name = T-Model-LandLab
68+
set Script argument =
69+
set Script path =
70+
71+
set Hillslope diffusion coefficient = 1e-4
72+
set Stream power erodibility coefficient = 1e-5
73+
set Define LandLab parameters using years = true
74+
75+
# 2.5 km mesh for LandLab.
76+
subsection Mesh information
77+
set X extent = 102.5e3
78+
set Y extent = 100e3
79+
set Spacing = 2500
80+
end
81+
end
82+
end
83+
84+
# Vertical gravity
85+
subsection Gravity model
86+
set Model name = vertical
87+
88+
subsection Vertical
89+
set Magnitude = 0
90+
end
91+
end
92+
93+
# One material with unity properties
94+
subsection Material model
95+
set Model name = simple
96+
97+
subsection Simple model
98+
set Reference density = 1
99+
set Reference specific heat = 1
100+
set Reference temperature = 0
101+
set Thermal conductivity = 1
102+
set Thermal expansion coefficient = 0
103+
set Viscosity = 1
104+
end
105+
end
106+
107+
# We also have to specify that we want to use the Boussinesq
108+
# approximation (assuming the density in the temperature
109+
# equation to be constant, and incompressibility).
110+
subsection Formulation
111+
set Formulation = Boussinesq approximation
112+
end
113+
114+
# We here use a globally refined mesh without
115+
# adaptive mesh refinement.
116+
subsection Mesh refinement
117+
set Initial global refinement = 0
118+
set Initial adaptive refinement = 0
119+
set Time steps between mesh refinement = 0
120+
end
121+
122+
subsection Postprocess
123+
set List of postprocessors = visualization, topography, velocity statistics
124+
125+
subsection Topography
126+
set Output to file = true
127+
set Time between text output = 50e3
128+
end
129+
130+
subsection Visualization
131+
set Time between graphical output = 50e3
132+
set Output mesh velocity = true
133+
set Output mesh displacement = true
134+
set Output undeformed mesh = false
135+
set Interpolate output = false
136+
end
137+
138+
end

‎source/mesh_deformation/external_tool_interface.cc‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -42,8 +42,9 @@ namespace aspect
4242
AffineConstraints<double> &mesh_velocity_constraints,
4343
const std::set<types::boundary_id> &boundary_ids) const
4444
{
45-
// First compute a (global) vector that has the correct velocities
46-
// set at all boundary nodes:
45+
// First compute a (global) vector that has the correct velocities set at all boundary nodes.
46+
// This takes the points on the LandLab mesh (which are the evaluation points) and interpolates
47+
// the solution from ASPECT onto these evaluation points.
4748
const std::vector<std::vector<double>> aspect_surface_velocities = evaluate_aspect_variables_at_points();
4849

4950
const std::vector<Tensor<1,dim>> external_surface_velocities

0 commit comments

Comments
 (0)