Skip to content

Commit fb77803

Browse files
matthewhoffmanmatthewhoffman
authored
Merge pull request #728 from trhille/landice/convert_exo_to_mpas_molho
Add ability to convert exo to MPAS mesh when using MOLHO for optimization Treat MALI layers as repeated instances of the Albany layer when necessary, so that we can convert an exodus file to a MALI .nc file when using the single-layer MOLHO solver for the optimization.
2 parents 1968e01 + 6efb772 commit fb77803

1 file changed

Lines changed: 158 additions & 28 deletions

File tree

‎landice/mesh_tools_li/conversion_exodus_init_to_mpasli_mesh.py‎

Lines changed: 158 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -58,18 +58,27 @@
5858

5959
# Map and copy Exodus data to MPAS data
6060

61-
# Create dictionary of variables that are supported by the script
62-
mpas_exodus_var_dic = {"beta":"basal_friction_log",
63-
"muFriction":"mu_log", \
64-
"stiffnessFactor":"stiffening_factor", \
65-
"uReconstructX":"solution_1", \
66-
"uReconstructY":"solution_2", \
67-
"temperature":"temperature", \
68-
"basalTemperature":"temperature", \
69-
"surfaceTemperature":"temperature", \
70-
"thickness":"ice_thickness"}
71-
#"surfaceAirTemperature":"surface_air_temperature", \ #use with caution!
72-
# A mapping between mpas and exodus file. Add one if you need to manipulate a different variable!
61+
# Create dictionary mapping exodus variable names to MPAS variable names.
62+
# This allows us to check what fields are present in the exodus file and
63+
# map them to the correct MALI field name.
64+
exodus_to_mpas_dic = {"basal_friction_log":"beta",
65+
"mu_log":"muFriction",
66+
"stiffening_factor":"stiffnessFactor",
67+
"stiffening_factor_log":"stiffnessFactor",
68+
"solution_1":"uReconstructX",
69+
"solution_2":"uReconstructY",
70+
"temperature":"temperature",
71+
"ice_thickness":"thickness"}
72+
73+
# Build reverse lookup: MPAS variable name -> list of possible exodus names
74+
mpas_to_exodus_candidates = {}
75+
for exo_name, mpas_name in exodus_to_mpas_dic.items():
76+
if mpas_name not in mpas_to_exodus_candidates:
77+
mpas_to_exodus_candidates[mpas_name] = []
78+
mpas_to_exodus_candidates[mpas_name].append(exo_name)
79+
# surfaceTemperature and basalTemperature also use the "temperature" exodus field
80+
mpas_to_exodus_candidates["surfaceTemperature"] = ["temperature"]
81+
mpas_to_exodus_candidates["basalTemperature"] = ["temperature"]
7382
ice_density = 910.0
7483
ocean_density = 1028.0
7584

@@ -104,8 +113,18 @@
104113
# Albany layering is in reverse order from mpas
105114
exo_layer_thick = np.flip(np.reshape(np.array(exo_layer_thick), [len(exo_layer_thick),]))
106115
mpas_layer_thick = np.ma.filled(dataset.variables['layerThicknessFractions'][:])
107-
if mpas_layer_thick.any() != exo_layer_thick.any():
116+
nLayers_exo = len(exo_layer_thick)
117+
nLayers_mpas = len(mpas_layer_thick)
118+
119+
if nLayers_exo <= 1 and nLayers_mpas > 1:
120+
print("WARNING: Albany exodus file has {} layer(s) but MPAS file has {} layers.".format(nLayers_exo, nLayers_mpas))
121+
print(" Assuming MOLHO (mono-layer higher-order) Albany solve.")
122+
print(" 3D fields will be interpolated to MPAS vertical layers using sigma^4 profile.")
123+
molho_mode = True
124+
elif np.array_equal(mpas_layer_thick, exo_layer_thick) == False:
108125
sys.exit("ERROR: Albany layer_thickness_ratio does not match MPAS layerThicknessFractions! Aborting")
126+
else:
127+
molho_mode = False
109128

110129
# Read cellID_array from ascii file
111130
if not options.id_file:
@@ -161,22 +180,40 @@
161180
extrapolation = None
162181

163182
print("\n---------------")
164-
if not var_name in mpas_exodus_var_dic:
165-
sys.exit("ERROR: variable '{}' not supported by this script. Supported variables are: {}".format(var_name, mpas_exodus_var_dic.keys()))
166-
if mpas_exodus_var_dic[var_name] in exo.get_node_variable_names() and var_name in dataset.variables.keys():
167-
exo_var_name = mpas_exodus_var_dic[var_name]
183+
if not var_name in mpas_to_exodus_candidates:
184+
sys.exit("ERROR: variable '{}' not supported by this script. Supported variables are: {}".format(var_name, mpas_to_exodus_candidates.keys()))
185+
186+
# Find the exodus variable name corresponding to this MPAS variable
187+
exo_var_name = None
188+
for candidate in mpas_to_exodus_candidates[var_name]:
189+
if candidate in exo.get_node_variable_names():
190+
exo_var_name = candidate
191+
break
192+
193+
# In MOLHO mode, temperature is derived from enthalpy (solution_4) in the exo file
194+
if molho_mode and var_name == "temperature":
195+
if "solution_4" in exo.get_node_variable_names():
196+
exo_var_name = "solution_4"
197+
print("MOLHO mode: reading enthalpy (solution_4) from exo to derive temperature")
198+
else:
199+
sys.exit("ERROR: MOLHO mode requires 'solution_4' (enthalpy) in the exo file for temperature conversion.")
200+
201+
if exo_var_name is not None and var_name in dataset.variables.keys():
168202
data_exo = np.array(exo.get_node_variable_values(exo_var_name,1))
169203

170204
# Get number of vertical layers from mpas output file.
171205
if len(np.shape(dataset.variables[var_name])) == 3:
172-
if var_name == "temperature":
206+
if molho_mode:
207+
nVert_max = 2 # MOLHO exo has 2 interfaces (top and bottom) for 3D fields
208+
elif var_name == "temperature":
173209
nVert_max = np.shape(dataset.variables[var_name])[2] + 1 #albany temperature is on diferent vertical grid than MPAS, and has one extra layer
174210
else:
175211
nVert_max = np.shape(dataset.variables[var_name])[2]
176212
else:
177213
nVert_max = 1
178214

179215
albanyTemperature = np.zeros((nCells, nVert_max)) #initialize albanyTemperature array to fill in below
216+
albanyVelocity = np.zeros((nCells, nVert_max)) #initialize albanyVelocity array for MOLHO sigma^4 interpolation
180217
dataset.variables[var_name][:] = 0 # Fill variable with zeros in order to ensure proper values in void
181218

182219
print("Begin {} conversion".format(var_name))
@@ -185,7 +222,9 @@
185222
print('Converting layer/level {} of {}'.format(nVert+1, nVert_max))
186223

187224
#Albany has inverted layer/level ordering relative to MPAS.
188-
if var_name == "surfaceTemperature":
225+
if molho_mode and var_name == "surfaceTemperature":
226+
nVert_albany = 1 # Top interface of MOLHO exo mesh (2 interfaces: 0=bed, 1=top)
227+
elif var_name == "surfaceTemperature":
189228
nVert_albany = nVert_max - 1
190229
elif var_name == "basalTemperature":
191230
nVert_albany = 0 # This not needed, because it will be 0 for basalTemperature, but adding this case to make that assumption explicit
@@ -220,7 +259,10 @@
220259
elif var_name == "stiffnessFactor":
221260
dataset.variables[var_name][0,cellID_array-1] = np.exp(data_exo_layer)
222261
elif var_name == "uReconstructX" or var_name == "uReconstructY":
223-
dataset.variables[var_name][0,cellID_array-1, nVert] = data_exo_layer / (60. * 60. * 24 * 365)
262+
if molho_mode:
263+
albanyVelocity[cellID_array-1, nVert] = data_exo_layer / (60. * 60. * 24 * 365)
264+
else:
265+
dataset.variables[var_name][0,cellID_array-1, nVert] = data_exo_layer / (60. * 60. * 24 * 365)
224266
elif var_name == "thickness":
225267
print("WARNING: thickness conversion is still experimental! Carefully check results before using.")
226268
# We have to be careful to not change MPAS geometry in the extended layer of cells - only touch it where the MPAS file already had nonzero thickness
@@ -242,14 +284,85 @@
242284
dataset.variables[var_name][0,cellID_array-1] = data_exo_layer
243285

244286
if var_name == "temperature":
245-
# Make generic equally-spaced layers for albany and MPAS. This works
246-
# even for unequally-spaced layers because this is a linear interpolation
247-
# and MPAS temperature layers are always halfway between Albany temperature
248-
# layers.
249-
albany_layers = np.arange(0, nVert_max)
250-
MPAS_layers = np.arange(1, nVert_max) - 0.5 # MPAS and albany temperature layers are staggered
251-
temperatureInterpolant = interp1d(albany_layers, albanyTemperature, axis=1)
252-
dataset.variables["temperature"][0,:,:] = temperatureInterpolant(MPAS_layers)
287+
if molho_mode:
288+
# In MOLHO mode, Albany has 2 temperature interfaces (top and bottom).
289+
# The sigma^4 profile applies to ENTHALPY, not temperature directly:
290+
# E(sigma) = E_bed + (E_top - E_bed) * (1 - sigma^4)
291+
# T_melt(sigma) = Tm - CCcoeff * ice_density * g * H * sigma
292+
# T(sigma) = min(1e6 / (ice_density * c_i) * E(sigma) + T0, T_melt(sigma))
293+
# where sigma is 0 at the top surface and 1 at the bed.
294+
# albanyTemperature[:, 0] = T_top (extracted from Albany top interface)
295+
# albanyTemperature[:, 1] = T_bed (extracted from Albany bottom interface)
296+
297+
# Physical constants (should match Albany input YAML "LandIce Physical Parameters")
298+
c_i = 2009.0 # "Heat capacity of ice" [J/(kg*K)]
299+
g = 9.80616 # "Gravity Acceleration" [m/s^2]
300+
T0 = 265.0 # "Reference Temperature" [K]
301+
CCcoeff = 7.9e-08 # "Clausius-Clapeyron Coefficient" [K/Pa]
302+
Tm = 273.15 # "Atmospheric Pressure Melting Temperature" [K]
303+
304+
# albanyTemperature here actually contains enthalpy values
305+
# read directly from "solution_4" in the exo file
306+
E_top = albanyTemperature[:, 0]
307+
E_bed = albanyTemperature[:, 1]
308+
309+
# Get ice thickness in meters (from MPAS file)
310+
H = dataset.variables['thickness'][0, :]
311+
312+
# Verify enthalpy-to-temperature conversion against Albany's temperature field
313+
if "temperature" in exo.get_node_variable_names():
314+
data_temp_exo = np.array(exo.get_node_variable_values("temperature", 1))
315+
if ordering == 1.0:
316+
T_top_exo = data_temp_exo[1::int(stride)] # Albany top interface
317+
T_bed_exo = data_temp_exo[0::int(stride)] # Albany bottom interface
318+
elif ordering == 0.0:
319+
node_num_check = int(stride)
320+
T_bed_exo = data_temp_exo[0:node_num_check]
321+
T_top_exo = data_temp_exo[node_num_check:2*node_num_check]
322+
# Use Albany's thickness for pressure melting calculation (exo stores in km)
323+
if "ice_thickness" in exo.get_node_variable_names():
324+
data_thk_exo = np.array(exo.get_node_variable_values("ice_thickness", 1))
325+
if ordering == 1.0:
326+
H_albany = data_thk_exo[0::int(stride)] * 1000.0 # convert km to m; thickness same at all levels
327+
elif ordering == 0.0:
328+
H_albany = data_thk_exo[0:node_num_check] * 1000.0
329+
else:
330+
print(" WARNING: 'ice_thickness' not found in exo; using MPAS thickness for verification.")
331+
H_albany = H[cellID_array-1]
332+
# Compute T from E at top (sigma=0) and bed (sigma=1)
333+
T_from_E_top = np.minimum(1.0e6 / (ice_density * c_i) * E_top[cellID_array-1] + T0, Tm)
334+
T_melt_bed = Tm - CCcoeff * ice_density * g * H_albany * 1.0
335+
T_from_E_bed = np.minimum(1.0e6 / (ice_density * c_i) * E_bed[cellID_array-1] + T0, T_melt_bed)
336+
max_diff_top = np.max(np.abs(T_from_E_top - T_top_exo))
337+
max_diff_bed = np.max(np.abs(T_from_E_bed - T_bed_exo))
338+
print(" Verification: max |T_from_enthalpy - T_exo| at top = {:.6e} K".format(max_diff_top))
339+
print(" Verification: max |T_from_enthalpy - T_exo| at bed = {:.6e} K".format(max_diff_bed))
340+
if max_diff_top > 0.01 or max_diff_bed > 0.01:
341+
print(" WARNING: Enthalpy-derived temperature differs from Albany temperature by > 0.01 K.")
342+
print(" Check that physical constants match the Albany input YAML.")
343+
344+
nVertLevels = dataset.dimensions["nVertLevels"].size
345+
layerThicknessFractions = np.ma.filled(dataset.variables['layerThicknessFractions'][:])
346+
print("MOLHO mode: interpolating temperature via enthalpy sigma^4 profile")
347+
for nLayer in np.arange(0, nVertLevels):
348+
# sigma at layer midpoint
349+
sigma = np.sum(layerThicknessFractions[0:nLayer]) + layerThicknessFractions[nLayer] / 2.0
350+
# Interpolate enthalpy
351+
E_sigma = E_bed[cellID_array-1] + (E_top[cellID_array-1] - E_bed[cellID_array-1]) * (1.0 - sigma**4)
352+
# Pressure melting temperature at this depth
353+
T_melt = Tm - CCcoeff * ice_density * g * H[cellID_array-1] * sigma
354+
# Convert enthalpy to temperature, capped at pressure melting point
355+
dataset.variables["temperature"][0,cellID_array-1,nLayer] = np.minimum(
356+
1.0e6 / (ice_density * c_i) * E_sigma + T0, T_melt)
357+
else:
358+
# Make generic equally-spaced layers for albany and MPAS. This works
359+
# even for unequally-spaced layers because this is a linear interpolation
360+
# and MPAS temperature layers are always halfway between Albany temperature
361+
# layers.
362+
albany_layers = np.arange(0, nVert_max)
363+
MPAS_layers = np.arange(1, nVert_max) - 0.5 # MPAS and albany temperature layers are staggered
364+
temperatureInterpolant = interp1d(albany_layers, albanyTemperature, axis=1)
365+
dataset.variables["temperature"][0,:,:] = temperatureInterpolant(MPAS_layers)
253366

254367
# Add reasonable (non-zero) temperatures outside of ice mask. Make this a linear
255368
# interpolation between min(273.15, surfaceAirTemperature) at the surface and 268.15K at the bed.
@@ -263,6 +376,23 @@
263376

264377
print('\nTemperature interpolation complete')
265378

379+
# In MOLHO mode, interpolate velocity to all MPAS levels using sigma^4 profile
380+
if molho_mode and var_name in ("uReconstructX", "uReconstructY"):
381+
# albanyVelocity[:, 0] = u_top (extracted from Albany top interface)
382+
# albanyVelocity[:, 1] = u_bed (extracted from Albany bottom interface)
383+
# u(sigma) = u_bed + (1 - sigma^4) * (u_top - u_bed)
384+
# where sigma is 0 at the top surface and 1 at the bed.
385+
u_top = albanyVelocity[:, 0]
386+
u_bed = albanyVelocity[:, 1]
387+
nVertLevels = dataset.dimensions["nVertLevels"].size
388+
nVertInterfaces = nVertLevels + 1
389+
layerThicknessFractions = np.ma.filled(dataset.variables['layerThicknessFractions'][:])
390+
print("MOLHO mode: interpolating {} using sigma^4 profile".format(var_name))
391+
for iInterface in np.arange(0, nVertInterfaces):
392+
# sigma at interface position (0 at top, 1 at bed)
393+
sigma = np.sum(layerThicknessFractions[0:iInterface])
394+
dataset.variables[var_name][0,:,iInterface] = u_bed + (1.0 - sigma**4) * (u_top - u_bed)
395+
266396
# Extrapolate and smooth beta and stiffnessFactor fields
267397
if var_name in ["beta", "muFriction", "stiffnessFactor"]:
268398
# Extrapolation

0 commit comments

Comments
 (0)