Skip to content

Commit 003a072

Browse files
danieldouglas92danieldouglas92
authored
Merge pull request #27 from tjhei/lla-spherical
spherical landlab
2 parents 3df4a55 + 9657355 commit 003a072

File tree

6 files changed

+365
-19
lines changed

6 files changed

+365
-19
lines changed

‎cookbooks/landlab/spherical/spherical_landlab.py‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
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 cartesian_to_spherical(x, y, z):
28+
radius = np.sqrt(x**2 + y**2 + z**2)
29+
theta = np.rad2deg(np.arccos(z / radius))
30+
phi = np.rad2deg(np.arctan2(y, x))
31+
32+
return radius, theta, phi
33+
34+
def initialize(comm_handle):
35+
if not comm_handle is None:
36+
# Convert the handle back to an MPI communicator
37+
global comm
38+
comm = MPI.Comm.f2py(comm_handle)
39+
40+
rank = comm.Get_rank()
41+
size = comm.Get_size()
42+
43+
print(f"Python: Hello from Rank {rank} of {size}")
44+
45+
data = 1
46+
globalsum = comm.allreduce(data, op=MPI.SUM)
47+
if comm.rank == 0:
48+
print(f"\tPython: testing communication; sum {globalsum}")
49+
else:
50+
print("Python: running sequentially!")
51+
52+
def finalize():
53+
pass
54+
55+
56+
# Run the Landlab simulation from the current time to end_time and return
57+
# the new topographic elevation (in m) at each local node.
58+
# dict_variable_name_to_value_in_nodes is a dictionary mapping variables
59+
# (x velocity, y velocity, temperature, etc.) to an array of values in each
60+
# node.
61+
def update_until(end_time, dict_variable_name_to_value_in_nodes):
62+
global current_time, elevation, linear_diffuser, flow_accumulator, stream_power_eroder, timestep
63+
64+
dt = end_time - current_time
65+
timestep += 1
66+
67+
deposition_erosion = np.zeros(model_grid.number_of_nodes)
68+
69+
x_velocity = dict_variable_name_to_value_in_nodes["x velocity"]
70+
y_velocity = dict_variable_name_to_value_in_nodes["y velocity"]
71+
z_velocity = dict_variable_name_to_value_in_nodes["z velocity"]
72+
73+
radial_velocity = np.sqrt(x_velocity**2 + y_velocity**2 + z_velocity**2)
74+
75+
# Substepping for surface processes
76+
if dt>0:
77+
n_substeps = 10
78+
sub_dt = dt / n_substeps
79+
for _ in range(n_substeps):
80+
81+
# TODO:
82+
elevation_before = elevation.copy()
83+
84+
#flow_accumulator.run_one_step()
85+
#stream_power_eroder.run_one_step(sub_dt)
86+
linear_diffuser.run_one_step(sub_dt)
87+
88+
elevation += radial_velocity * sub_dt
89+
90+
deposition_erosion += elevation - elevation_before
91+
pass
92+
93+
filename = f"./output-spherical-landlab/landlab_{str(timestep).zfill(3)}.vtk"
94+
print("Writing output VTK file...", filename)
95+
vtk_file = write_legacy_vtk(path=filename, grid=model_grid, clobber=True, z_at_node=elevation)
96+
97+
vtks.append((current_time, filename))
98+
99+
if True:
100+
# write vtk.series file (ParaView supports legacy VTK in this format)
101+
with open("./output-spherical-landlab/landlab.vtk.series", "w") as f:
102+
series = {
103+
"file-series-version": "1.0",
104+
"files": [
105+
{"name": os.path.basename(filename), "time": time / s2yr}
106+
for time, filename in vtks
107+
]
108+
}
109+
json.dump(series, f, indent=2)
110+
111+
current_time = end_time
112+
print("Max elevation:", np.max(elevation),
113+
"Min elevation:", np.min(elevation),
114+
"min deposition:", np.min(deposition_erosion),
115+
"max deposition:", np.max(deposition_erosion))
116+
117+
return deposition_erosion.copy()
118+
119+
def set_mesh_information(dict_grid_information):
120+
global model_grid, elevation
121+
122+
if not model_grid:
123+
print("* Creating Spherical Grid ...")
124+
#x_extent = dict_grid_information["Mesh X extent"]
125+
#y_extent = dict_grid_information["Mesh Y extent"]
126+
#spacing = dict_grid_information["Mesh Spacing"]
127+
128+
model_grid = landlab.IcosphereGlobalGrid(radius=6371e3*0.99, mesh_densification_level=3)
129+
print(type(model_grid.x_of_node),type(model_grid.y_of_node),type(model_grid.z_of_node))
130+
for i in range(len(model_grid.x_of_node)):
131+
x=model_grid.x_of_node[i]
132+
y=model_grid.y_of_node[i]
133+
z=model_grid.z_of_node[i]
134+
print(i,x,y,z,np.sqrt(x*x+y*y+z*z))
135+
136+
137+
138+
139+
print("model grid nodes", len(model_grid.x_of_node))
140+
r, theta, phi = cartesian_to_spherical(model_grid.x_of_node, \
141+
model_grid.y_of_node, \
142+
model_grid.z_of_node)
143+
144+
elevation = model_grid.add_zeros("topographic__elevation", at="node")
145+
elevation += 1e3 + np.random.rand(elevation.size) * 0.1e3
146+
#elevation[theta < 90] += 1000
147+
print(elevation)
148+
149+
initialize_landlab_components(None)
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_dictionary):
155+
global model_grid
156+
#dim = grid_dictionary["ASPECT Dimension"]
157+
#if dim == 2:
158+
# print(np.unique(model_grid.x_of_node))
159+
# return np.unique(model_grid.x_of_node)
160+
#if dim == 3:
161+
return model_grid.x_of_node.copy()
162+
163+
# Return the y coordinates of the locally owned nodes on this
164+
# MPI rank. grid_id is always 0.
165+
def get_grid_y(grid_dictionary):
166+
global model_grid
167+
return model_grid.y_of_node.copy()
168+
def get_grid_z(grid_dictionary):
169+
global model_grid
170+
return model_grid.z_of_node.copy()
171+
172+
def initialize_landlab_components(landlab_component_parameters):
173+
global model_grid, elevation, linear_diffuser, flow_accumulator, stream_power_eroder
174+
175+
D = 1e-4 #landlab_component_parameters["Hillslope coefficient"]
176+
K_sp = 1e-5 #landlab_component_parameters["Stream power erodibility coefficient"]
177+
flow_director = "Steepest" # Only flow director supported for HexGrids
178+
179+
180+
print("* Creating LinearDiffuser ... with D =", D)
181+
linear_diffuser = LinearDiffuser(model_grid, linear_diffusivity=D)
182+
print("* Creating FlowAccumulator ... with flow director =", flow_director)
183+
flow_accumulator = FlowAccumulator(model_grid, elevation, flow_director=flow_director)
184+
print("* Creating StreamPowerEroder ... with K_sp =", K_sp)
185+
stream_power_eroder = StreamPowerEroder(model_grid, K_sp=K_sp)
186+
187+
pass
188+
# Return the initial topography at the start of the simulation
189+
# in each node.
190+
def get_initial_topography(grid_dictionary):
191+
global elevation
192+
#dim = grid_dictionary["ASPECT Dimension"]
193+
#if dim == 2:
194+
# return elevation[model_grid.y_of_node == 0]
195+
#if dim == 3:
196+
return elevation
197+

‎cookbooks/landlab/spherical/spherical_test.prm‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
# This parameter file sets up a model to test the interfacing between
2+
# the surface processes code LandLab and ASPECT. The model prescribes
3+
# uniform uplift within the ASPECT model (1 cm/yr), which is then used
4+
# to uplift the LandLab model. The LandLab model erodes the topography
5+
# using Hillslope Diffusion and Stream power.
6+
set Dimension = 3
7+
8+
9+
set Use years instead of seconds = true
10+
set End time = 2000e3
11+
set Maximum time step = 10e1
12+
13+
14+
#set Nonlinear solver scheme = single Advection, no Stokes
15+
set Nonlinear solver scheme = no Advection, no Stokes
16+
set Pressure normalization = surface
17+
set Surface pressure = 0
18+
set Output directory = output-spherical-landlab
19+
20+
subsection Geometry model
21+
set Model name = spherical shell
22+
23+
subsection Spherical shell
24+
25+
set Inner radius = 3481e3
26+
set Outer radius = 6371e3
27+
28+
end
29+
end
30+
31+
# Temperature effects are ignored
32+
subsection Initial temperature model
33+
set Model name = function
34+
subsection Function
35+
set Function expression = 273
36+
end
37+
end
38+
39+
subsection Boundary temperature model
40+
set Fixed temperature boundary indicators = bottom, top
41+
set List of model names = initial temperature
42+
end
43+
44+
# Free slip on all boundaries
45+
subsection Boundary velocity model
46+
set Tangential velocity boundary indicators = bottom
47+
end
48+
49+
subsection Mesh deformation
50+
set Mesh deformation boundary indicators = top: Landlab
51+
52+
subsection Landlab
53+
set MPI ranks for Landlab = 1
54+
set Script name = spherical_landlab
55+
set Script argument =
56+
set Script path =
57+
end
58+
end
59+
60+
# Vertical gravity
61+
subsection Gravity model
62+
set Model name = radial constant
63+
64+
subsection Radial constant
65+
set Magnitude = 1
66+
end
67+
end
68+
69+
# One material with unity properties
70+
subsection Material model
71+
set Model name = simple
72+
73+
subsection Simple model
74+
set Reference density = 1
75+
set Reference specific heat = 1
76+
set Reference temperature = 0
77+
set Thermal conductivity = 1
78+
set Thermal expansion coefficient = 0
79+
set Viscosity = 1
80+
end
81+
end
82+
83+
# We also have to specify that we want to use the Boussinesq
84+
# approximation (assuming the density in the temperature
85+
# equation to be constant, and incompressibility).
86+
subsection Formulation
87+
set Formulation = Boussinesq approximation
88+
end
89+
90+
# We here use a globally refined mesh without
91+
# adaptive mesh refinement.
92+
subsection Mesh refinement
93+
set Initial global refinement = 0
94+
set Initial adaptive refinement = 0
95+
set Time steps between mesh refinement = 0
96+
end
97+
98+
subsection Postprocess
99+
set List of postprocessors = visualization, topography, velocity statistics
100+
101+
subsection Topography
102+
set Output to file = true
103+
set Time between text output = 50e3
104+
end
105+
106+
subsection Visualization
107+
set Time between graphical output = 50e3
108+
set Output mesh velocity = true
109+
set Output mesh displacement = true
110+
set Output undeformed mesh = false
111+
set Interpolate output = false
112+
end
113+
114+
end

‎include/aspect/mesh_deformation/landlab.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -121,6 +121,10 @@ namespace aspect
121121
*/
122122
std::string script_argument;
123123

124+
/**
125+
* Whether the ASPECT geometry is spherical.
126+
*/
127+
bool is_spherical;
124128
#ifdef ASPECT_WITH_PYTHON
125129
/**
126130
* The Python module object.

‎source/mesh_deformation/external_tool_interface.cc‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@
2424

2525
#include <deal.II/numerics/vector_tools_evaluate.h>
2626

27+
2728
/*
2829
TODO:
2930
- initial topography from external tool: for now compute_initial_deformation_on_boundary(), later refactor
@@ -104,7 +105,9 @@ namespace aspect
104105
// so we need to use a simple mapping instead of the one stored in the Simulator. The latter
105106
// would produce the deformed mesh. We currently always use a Q1 mapping when mesh deformation
106107
// is enabled, so a Q1 mapping is the right choice.
107-
static MappingQ<dim> mapping(1);
108+
// TODO: if the mesh is cartesian, we could just use a Q1 mapping, but we need a higher order
109+
// mapping for spherical meshes. Which order should we pick?
110+
static MappingQ<dim> mapping(4);
108111
remote_point_evaluator = std::make_unique<Utilities::MPI::RemotePointEvaluation<dim, dim>>();
109112
remote_point_evaluator->reinit(this->evaluation_points, this->get_triangulation(), mapping);
110113

‎source/mesh_deformation/interface.cc‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1264,7 +1264,7 @@ namespace aspect
12641264
= sim.parameters.linear_stokes_solver_tolerance
12651265
* ((this->simulator_is_past_initialization()) ? 1.0 : 1e-5);
12661266

1267-
SolverControl solver_control_mf(5 * rhs.size(),
1267+
SolverControl solver_control_mf(200,
12681268
tolerance * rhs.l2_norm());
12691269
SolverCG<dealii::LinearAlgebra::distributed::Vector<double>> cg(solver_control_mf);
12701270

0 commit comments

Comments
 (0)