Skip to content

update for ALEIB_CM #754

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
julesghub merged 1 commit into from
May 12, 2026
Merged

update for ALEIB_CM #754

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
197 changes: 197 additions & 0 deletions docs/UWGeodynamics/examples/1_23_08_CouplingModellingALEIB_TopographyRelaxation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,197 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "9182138b-b954-4a5c-a6c6-35cb6c296923",
"metadata": {},
"source": [
"Topography Relaxation with surface processes\n",
"======\n",
"\n",
"This notebook models the topography relaxation with surface processes (Model TR-CM2) as described in Neng et al. (2026). It integrates Underworld2 with Badlands using the UWGeodynamic module within the ALE-IB scheme.\n",
"\n",
"**References**\n",
"\n",
"Lu, N., Moresi, L., Giordani, J., & Knight, B. (2026). A novel ALE scheme with the internal boundary for coupling tectonic and surface processes in geodynamic models. EGUsphere, 2026, 1-35."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "72d4536d-33cb-4c16-8812-645165b311f8",
"metadata": {},
"outputs": [],
"source": [
"import underworld as uw\n",
"import math\n",
"from underworld import function as fn\n",
"import numpy as np\n",
"import os\n",
"\n",
"from underworld import UWGeodynamics as GEO\n",
"u = GEO.UnitRegistry\n",
"ndim = GEO.non_dimensionalise\n",
"dimen = GEO.dimensionalise\n",
"\n",
"comm = uw.mpi.comm\n",
"rank = uw.mpi.rank\n",
"size = uw.mpi.size\n",
"\n",
"GEO.rcParams[\"initial.nonlinear.tolerance\"] = 1e-3\n",
"GEO.rcParams['initial.nonlinear.max.iterations'] = 100\n",
"GEO.rcParams[\"nonlinear.tolerance\"] = 1e-3\n",
"GEO.rcParams['nonlinear.max.iterations'] = 100\n",
"GEO.rcParams[\"popcontrol.particles.per.cell.2D\"] = 20\n",
"GEO.rcParams[\"swarm.particles.per.cell.2D\"] = 20\n",
"GEO.rcParams[\"surface.pressure.normalization\"] = True\n",
"GEO.rcParams[\"pressure.smoothing\"] = True\n",
"GEO.rcParams[\"popcontrol.split.threshold\"] = 0.1\n",
"\n",
"half_rate = 1.0 * u.centimeter / u.year\n",
"model_length = 500. * u.kilometer\n",
"gravity = 9.81 * u.meter / u.second**2\n",
"bodyforce = 3300 * u.kilogram / u.metre**3 *gravity \n",
"\n",
"KL = model_length\n",
"Kt = KL / half_rate\n",
"KM = bodyforce * KL**2 * Kt**2\n",
"\n",
"GEO.scaling_coefficients[\"[length]\"] = KL\n",
"GEO.scaling_coefficients[\"[time]\"] = Kt\n",
"GEO.scaling_coefficients[\"[mass]\"]= KM"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d497849b-031f-4eee-af78-a5579245b274",
"metadata": {},
"outputs": [],
"source": [
"longtest = False\n",
"dy = ndim(10.0 * u.kilometer)\n",
"max_time = 11.0*u.kiloyear\n",
"\n",
"if \"UW_LONGTEST\" in os.environ or longtest:\n",
" dy = ndim(2.5 * u.kilometer)\n",
" max_time = 301.0*u.kiloyear\n",
"\n",
"xmin, xmax = ndim(-160 * u.kilometer), ndim(160 * u.kilometer)\n",
"ymin, ymax = ndim(-160 * u.kilometer), ndim(40 * u.kilometer)\n",
"yint = 0.\n",
" \n",
"dy = ndim(2.5 * u.kilometer)\n",
"dx = dy\n",
"xRes,yRes = int(np.around((xmax-xmin)/dx)),int(np.around((ymax-ymin)/dy))\n",
"yResa,yResb =int(np.around((ymax-yint)/dy)),int(np.around((yint-ymin)/dy))\n",
"\n",
"Model = GEO.Model(elementRes=(xRes, yRes),\n",
" minCoord=(xmin,ymin),\n",
" maxCoord=(xmax, ymax),\n",
" gravity=(0.0, -gravity),\n",
" periodic=(True, False))\n",
"Model.outputDir= \"1_23_08_CouplingModellingALEIB_TopographyRelaxation_yres{:n}\".format(yRes)\n",
"Model.minStrainRate = 1e-18 / u.second\n",
"\n",
"wRatio = 1\n",
"D = np.abs(ymin)\n",
"Lambda = D/wRatio\n",
"k = 2.0 * np.pi / Lambda\n",
"mu0 = ndim(1e21 * u.pascal * u.second)\n",
"g = ndim(gravity)\n",
"rho0 = ndim(3300* u.kilogram / u.metre**3)\n",
"drho = rho0-0.\n",
"w_m = ndim(5*u.kilometer)\n",
"\n",
"tau0 = 2*k*mu0/drho/g\n",
"tau = (D*k+np.sinh(D*k)*np.cosh(D*k))/(np.sinh(D*k)**2)*tau0\n",
"\n",
"fn_coord = fn.input()\n",
"surf_fn = w_m * fn.math.cos(2.*np.pi*fn_coord[0]/Lambda) \n",
"\n",
"Model.inter_wall = Model._get_InternalwallSets(yint)\n",
"with Model.mesh.deform_mesh():\n",
" Model.mesh.data[Model.inter_wall.data, 1] = surf_fn.evaluate(Model.inter_wall)[:,0]\n",
"Model._freeSurface_ALEIB = True \n",
"Model.freeSurface = True \n",
"Model._freeSurface.solve(0.)\n",
"\n",
"# recreate swarm as mesh deformed\n",
"import underworld as uw\n",
"from collections import OrderedDict \n",
"Model.swarm_variables = OrderedDict()\n",
"Model.swarm = uw.swarm.Swarm(mesh=Model.mesh, particleEscape=True)\n",
"Model.swarm.allow_parallel_nn = True \n",
"particlesPerCell = GEO.rcParams[\"swarm.particles.per.cell.2D\"]\n",
"Model._swarmLayout = uw.swarm.layouts.PerCellSpaceFillerLayout(swarm=Model.swarm,particlesPerCell=particlesPerCell)\n",
"Model.swarm.populate_using_layout(layout=Model._swarmLayout)\n",
"Model._initialize()\n",
"\n",
"materialAShape = fn_coord[1] > surf_fn\n",
"materialMShape = fn_coord[1] <= surf_fn\n",
"\n",
"materialA = Model.add_material(name=\"Air\", shape=materialAShape)\n",
"materialM = Model.add_material(name=\"Mantle\", shape=materialMShape) \n",
"\n",
"materialA.viscosity = 1e18 * u.pascal * u.second\n",
"materialM.viscosity = 1e21 * u.pascal * u.second\n",
"\n",
"materialA.density = 0.\n",
"materialM.density = 3300 * u.kilogram / u.metre**3\n",
"\n",
"sediment = Model.add_material(name=\"sediment\")\n",
"sediment.viscosity = 1e19 * u.pascal * u.second\n",
"sediment.density = 2700 * u.kilogram / u.metre**3 \n",
"\n",
"Model.set_velocityBCs(left=[0.,None],right=[0,None],bottom=[0.,0.], top=[None, 0.])\n",
"Model.init_model()\n",
"\n",
"Model.solver.set_inner_method(\"mumps\")\n",
"Model.solver.set_penalty(1e3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "72102da9-434f-4de3-b373-4c8db080694a",
"metadata": {},
"outputs": [],
"source": [
"dt_set = 2.5*u.kiloyear\n",
"checkpoint_interval = 10.0*u.kiloyear\n",
"Model.surfaceProcesses = GEO.surfaceProcesses.Badlands(airIndex=[materialA.index],sedimentIndex=sediment.index,XML=\"resources/badlands_1en5.xml\", resolution=1.25 * u.kilometre, checkpoint_interval=dt_set,aspectRatio2d=0.25,surfElevation=surf_fn)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "daf7ba37-4acc-4054-9843-f9526ac32a9a",
"metadata": {},
"outputs": [],
"source": [
"Model.run_for(max_time, checkpoint_interval=checkpoint_interval,dt=dt_set)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
111 changes: 111 additions & 0 deletions docs/UWGeodynamics/examples/resources/badlands_1en5.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
<?xml version="1.0" encoding="UTF-8"?>
<badlands xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<!-- Regular grid structure -->
<grid>
<!-- Optional parameter (integer) used to decrease TIN resolution.
The default value is set to 1. Increasing the factor
value will multiply the digital elevation model resolution
accordingly. -->
<resfactor>1</resfactor>
<boundary>flat</boundary>
<!-- <boundary>fixed</boundary> -->
<!-- Underworld flag. -->
<udw>1</udw>
</grid>
<!-- Sea-level structure -->
<sea>
<!-- Relative sea-level position [m] -->
<position>-0.</position>
<!-- Sea-level curve - (optional) -->
<!-- <curve>data/sea.csv</curve> -->
<!-- Limit flow network computation based on
water depth [m] -->
<!-- <limit>100.</limit> -->
</sea>
<!-- Simulation time structure -->
<time>
<!-- Simulation start time [a] -->
<start>0.</start>
<!-- Simulation end time [a] -->
<end>510000.</end>
<!-- Display interval [a] -->
<!-- This will be overridden by the linkage, so it doesn't matter -->
<display>10000.</display>
</time>
<strata>
<stratdx>1250.</stratdx>
<laytime>1250.</laytime>
</strata>
<!-- Precipitation structure -->
<precipitation>
<!-- Number of precipitation events -->
<climates>1</climates>
<!-- Precipitation definition -->
<rain>
<!-- Rain start time [a] -->
<rstart>0.</rstart>
<!-- Rain end time [a] -->
<rend>510000.</rend>
<!-- Rain computation time step [a] -->
<ortime>1250.</ortime>
<!-- Background precipitation value [m/a] -->
<rbgd>1.</rbgd>
<!-- Minimal precipitation value [m/a] -->
<rmin>0.1</rmin>
<!-- Maximal precipitation value [m/a] -->
<rmax>2.</rmax>
<!-- Wind velocity along X (W-E) direction [m/s] -->
<windx>2.</windx>
<!-- Wind velocity along Y (S-N) direction [m/s] -->
<windy>0.</windy>
<!-- Time conversion from cloud water to hydrometeors
range from 200 to 2000 [s]. Optional default is set
to 1000 s -->
<tauc>1000.</tauc>
<!-- Time for hydrometeor fallout range from 200 to 2000 [s].
Optional default is set to 1000 s -->
<tauf>1000.</tauf>
<!-- Moist stability frequency range from 0 to 0.01 [/s].
Optional default is set to 0.005 /s -->
<nm>0.005</nm>
<!-- Uplift sensitivity factor range from 0.001 to 0.02 [kg/m3].
Optional default is set to 0.005 kg/m3 -->
<cw>0.01</cw>
<!-- Depth of the moist layer range from 1000 to 5000 [m].
Optional default is set to 3000 m -->
<hw>5000.</hw>
</rain>
</precipitation>
<!-- Stream power law parameters:
The stream power law is a simplified form of the usual expression of
sediment transport by water flow, in which the transport rate is assumed
to be equal to the local carrying capacity, which is itself a function of
boundary shear stress. -->
<sp_law>
<!-- Values of m and n indicate how the incision rate scales
with bed shear stress for constant value of sediment flux
and sediment transport capacity.
Generally, m and n are both positive, and their ratio
(m/n) is considered to be close to 0.5 -->
<m>0.5</m>
<n>1.0</n>
<!-- The erodibility coefficient is scale-dependent and its value depend
on lithology and mean precipitation rate, channel width, flood
frequency, channel hydraulics. -->
<erodibility>1.e-5</erodibility>
<diffnb>10</diffnb>
<diffprop>0.05</diffprop>
</sp_law>
<!-- Linear slope diffusion parameters:
Parameterisation of the sediment transport includes the simple creep transport
law which states that transport rate depends linearly on topographic gradient. -->
<creep>
<!-- Surface diffusion coefficient [m2/a] -->
<caerial>0.1</caerial>
<!-- Marine diffusion coefficient [m2/a] -->
<cmarine>0.1</cmarine>
<criver>0.1</criver>
</creep>
<!-- Output folder path -->
<outfolder>outbdls</outfolder>
</badlands>
23 changes: 20 additions & 3 deletions src/underworld/UWGeodynamics/_freesurface.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,17 +100,34 @@ def _advect_surface(self, dt):
uw.mpi.barrier()
self.TField.syncronise()

def _update_mesh(self):
def _advect_surface_sp(self,dt):
surf_fn_badlands = self.model.surfaceProcesses.solve(dt)
if self.interface:
if self.model.mesh.dim == 2:
interpolate_x = self.model.mesh.data[self.interface.data,0]
interpolate_z = surf_fn_badlands(interpolate_x)
self.TField.data[self.interface.data, 0] = interpolate_z.copy()

if self.model.mesh.dim == 3:
interpolate_x = self.model.mesh.data[self.interface.data,0]
interpolate_x = self.model.mesh.data[self.interface.data,1]
interpolate_z = surf_fn_badlands((interpolate_x,interpolate_y))
self.TField.data[self.interface.data, 0] = interpolate_z.copy()
uw.mpi.barrier()
self.TField.syncronise()

def _update_mesh(self):
with self.model.mesh.deform_mesh():
# Last dimension is the vertical dimension
self.model.mesh.data[:, -1] = self.TField.data[:, 0].copy()

def solve(self, dtime):
""" Advect free surface through dt and update the mesh """

# First we advect the surface
self._advect_surface(dtime)
if self.model.surfaceProcesses:
self._advect_surface_sp(dtime)
else:
self._advect_surface(dtime)
# Then we solve the system of linear equation
self._solve_sle()
# Finally we update the mesh
Expand Down
Loading