Skip to content

Commit 809545b

Browse files
danieldouglas92danieldouglas92
committed
add example
1 parent 003a072 commit 809545b

File tree

2 files changed

+318
-0
lines changed

2 files changed

+318
-0
lines changed

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

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
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+
57+
dt = end_time - current_time
58+
timestep += 1
59+
60+
deposition_erosion = np.zeros(model_grid.number_of_nodes)
61+
62+
slice_x_velocity = dict_variable_name_to_value_in_nodes["x velocity"]
63+
slice_y_velocity = dict_variable_name_to_value_in_nodes["y 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+
vertical_velocity = np.zeros(model_grid.number_of_nodes)
74+
unique_x_values = np.unique(model_grid.x_of_node)
75+
for x in unique_x_values:
76+
vertical_velocity[model_grid.x_of_node == x] = slice_y_velocity[unique_x_values == x]
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 returning
118+
# the change in elevation used to deform the ASPECT mesh. Averaging is straight forward to
119+
# do going forward, and the commented line in the for loop below shows how this would be done.
120+
# However, for the purposes of this test I want to make sure the ASPECT mesh is
121+
# deforming to the exact same topography as the LandLab mesh where they intersect.
122+
deposition_erosion_2d = np.zeros(len(np.unique(model_grid.x_of_node)))
123+
for x in unique_x_values:
124+
# deposition_erosion_2d[unique_x_values == x] = np.average(deposition_erosion[model_grid.x_of_node == x])
125+
deposition_erosion_2d = deposition_erosion[model_grid.y_of_node == 0]
126+
127+
return deposition_erosion_2d
128+
129+
def set_mesh_information(dict_grid_information):
130+
global model_grid, elevation
131+
132+
if not model_grid:
133+
print("* Creating RasterModelGrid ...")
134+
x_extent = 100e3
135+
y_extent = 100e3
136+
spacing = 2500
137+
138+
nrows = int(y_extent / spacing) + 1 # number of node rows
139+
ncols = int(x_extent / spacing) + 1 # number of node columns
140+
141+
# For 2D models, shift the LandLab mesh so that it is centered on the y-axis.
142+
model_grid = landlab.RasterModelGrid((nrows, ncols), xy_spacing=(spacing, spacing), xy_of_lower_left=(0, -y_extent / 2))
143+
144+
model_grid.set_closed_boundaries_at_grid_edges(right_is_closed=True,
145+
left_is_closed=False,
146+
top_is_closed=True,
147+
bottom_is_closed=True)
148+
149+
print("* Creating topographic elevation ...")
150+
# Initialize topography array with zeros
151+
elevation = model_grid.add_zeros("topographic__elevation", at="node")
152+
# Assign 1000 m high topography with random noise
153+
np.random.seed(0) # For reproducibility
154+
elevation += 1000 + np.random.rand(model_grid.number_of_nodes)
155+
156+
initialize_landlab_components(None)
157+
print("* Done")
158+
159+
# Return the x coordinates of the locally owned nodes on this
160+
# MPI rank. grid_id is always 0.
161+
def get_grid_x(grid_dictionary):
162+
global model_grid
163+
return np.unique(model_grid.x_of_node)
164+
165+
# Return the y coordinates of the locally owned nodes on this
166+
# MPI rank. grid_id is always 0.
167+
def get_grid_y(grid_dictionary):
168+
global model_grid
169+
return np.zeros(np.unique(model_grid.x_of_node).shape)
170+
171+
def initialize_landlab_components(grid_dictionary):
172+
global model_grid, elevation, linear_diffuser, flow_accumulator, stream_power_eroder
173+
174+
D = 1e-4 / s2yr
175+
K_sp = 1e-5 / s2yr
176+
flow_director = "D8" # Only flow director supported for HexGrids
177+
178+
179+
print("* Creating LinearDiffuser ... with D =", D)
180+
linear_diffuser = LinearDiffuser(model_grid, linear_diffusivity=D)
181+
print("* Creating FlowAccumulator ... with flow director =", flow_director)
182+
flow_accumulator = FlowAccumulator(model_grid, elevation, flow_director=flow_director)
183+
print("* Creating StreamPowerEroder ... with K_sp =", K_sp)
184+
stream_power_eroder = StreamPowerEroder(model_grid, K_sp=K_sp)
185+
pass
186+
187+
# Return the initial topography at the start of the simulation
188+
# in each node.
189+
def get_initial_topography(grid_dictionary):
190+
global elevation
191+
return elevation[model_grid.y_of_node == 0]

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

Lines changed: 127 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,127 @@
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 resolution, 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+
end
71+
end
72+
73+
# Vertical gravity
74+
subsection Gravity model
75+
set Model name = vertical
76+
77+
subsection Vertical
78+
set Magnitude = 1e-100
79+
end
80+
end
81+
82+
# One material with unity properties
83+
subsection Material model
84+
set Model name = simple
85+
86+
subsection Simple model
87+
set Reference density = 1
88+
set Reference specific heat = 1
89+
set Reference temperature = 0
90+
set Thermal conductivity = 1
91+
set Thermal expansion coefficient = 0
92+
set Viscosity = 1
93+
end
94+
end
95+
96+
# We also have to specify that we want to use the Boussinesq
97+
# approximation (assuming the density in the temperature
98+
# equation to be constant, and incompressibility).
99+
subsection Formulation
100+
set Formulation = Boussinesq approximation
101+
end
102+
103+
# We here use a globally refined mesh without
104+
# adaptive mesh refinement.
105+
subsection Mesh refinement
106+
set Initial global refinement = 0
107+
set Initial adaptive refinement = 0
108+
set Time steps between mesh refinement = 0
109+
end
110+
111+
subsection Postprocess
112+
set List of postprocessors = visualization, topography, velocity statistics
113+
114+
subsection Topography
115+
set Output to file = true
116+
set Time between text output = 50e3
117+
end
118+
119+
subsection Visualization
120+
set Time between graphical output = 50e3
121+
set Output mesh velocity = true
122+
set Output mesh displacement = true
123+
set Output undeformed mesh = false
124+
set Interpolate output = false
125+
end
126+
127+
end

0 commit comments

Comments
 (0)