Skip to content

Commit

Permalink
Merge pull request #20 from equinor/timeseries
Browse files Browse the repository at this point in the history 
Timeseries
  • Loading branch information
@adamchengtkc
adamchengtkc committed Apr 19, 2024
2 parents b578172 + bc1f728 commit 93f8812
Show file tree
Hide file tree
Showing 10 changed files with 154 additions and 3,870 deletions.
2 changes: 1 addition & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ others/
__pycache__/
*.egg-info/
.eggs/

case/
.testmondata
.pytest_cache
.coverage
Expand Down
92 changes: 61 additions & 31 deletions warmth/mesh_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import dolfinx
from petsc4py import PETSc
import ufl
import sys
from scipy.interpolate import LinearNDInterpolator
from progress.bar import Bar
from warmth.build import single_node, Builder
Expand Down Expand Up @@ -254,7 +255,6 @@ def write_hexa_mesh_resqml( self, out_path, tti):
if (lid0>=-1) and (lid0<100):
hexa_to_keep.append(h)
lid_to_keep.append(lid0)
# cell_id_to_keep.append(self.node_index[i])
cell_id_to_keep.append(hex_data_nodeID[i])
minY = np.amin(np.array ( [x_original_order[hi,1] for hi in h] ))
if abs( self.node1D[hex_data_nodeID[i]].Y - minY)>1:
Expand Down Expand Up @@ -290,7 +290,8 @@ def write_hexa_mesh_resqml( self, out_path, tti):

from os import path
filename_hex = path.join(out_path, self.modelName+'_hexa_'+str(tti)+'.epc')
write_hexa_grid_with_properties(filename_hex, np.array(points_cached), hexa_renumbered, "hexamesh",
points_cached=np.array(points_cached)
write_hexa_grid_with_properties(filename_hex, points_cached, hexa_renumbered, "hexamesh",
np.array(T_per_vertex_keep), np.array(age_per_vertex_keep), poro0_per_cell, decay_per_cell, density_per_cell,
cond_per_cell, rhp_per_cell, lid_per_cell)
return filename_hex
Expand All @@ -305,11 +306,9 @@ def write_hexa_mesh_timeseries( self, out_path):
returns the filename (of the .epc file) that was written
"""
hexaHedra, hex_data_layerID, hex_data_nodeID = self.buildHexahedra(keep_padding=False)

hexa_to_keep = []
p_to_keep = set()
lid_to_keep = []
cond_per_cell = []
cell_id_to_keep = []
for i,h in enumerate(hexaHedra):
lid0 = hex_data_layerID[i]
Expand All @@ -320,39 +319,45 @@ def write_hexa_mesh_timeseries( self, out_path):
hexa_to_keep.append(h)
lid_to_keep.append(lid0)
cell_id_to_keep.append(hex_data_nodeID[i])
# k_cond_mean = []
for hi in h:
p_to_keep.add(hi)
# k_cond_mean.append(self.thermalCond.x.array[hi])
# cond_per_cell.append( np.mean(np.array(k_cond_mean)))

# poro0_per_cell = np.array( [ self.getSedimentPropForLayerID('phi', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ] )
# decay_per_cell = np.array( [ self.getSedimentPropForLayerID('decay', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
# density_per_cell = np.array( [ self.getSedimentPropForLayerID('solidus', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
# cond_per_cell = np.array( [ self.getSedimentPropForLayerID('k_cond', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
# rhp_per_cell = np.array( [ self.getSedimentPropForLayerID('rhp', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
# lid_per_cell = np.array(lid_to_keep)

points_cached_series = []
Temp_per_vertex_series = []

for tti in self.time_indices:
x_original_order, T_per_vertex = self.get_node_pos_and_temp(tti)
T_per_vertex_filt = [ T_per_vertex[i] for i in range(x_original_order.shape[0]) if i in p_to_keep ]
Temp_per_vertex_series.append(np.array(T_per_vertex_filt))
points_cached = []
point_original_to_cached = np.ones(x_original_order.shape[0], dtype = np.int32) * (-1)
for i in range(x_original_order.shape[0]):

poro0_per_cell = np.array( [ self.getSedimentPropForLayerID('phi', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ] )
decay_per_cell = np.array( [ self.getSedimentPropForLayerID('decay', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
density_per_cell = np.array( [ self.getSedimentPropForLayerID('solidus', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
cond_per_cell = np.array( [ self.getSedimentPropForLayerID('k_cond', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
rhp_per_cell = np.array( [ self.getSedimentPropForLayerID('rhp', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
lid_per_cell = np.array(lid_to_keep)

x_original_order, T_per_vertex = self.get_node_pos_and_temp(self.time_indices[0]) # oldest first
n_vertices = x_original_order.shape[0]
age_per_vertex_keep = np.array([ self.age_per_vertex[i] for i in range(n_vertices) if i in p_to_keep ])
Temp_per_vertex_series = np.empty([len(self.time_indices), len(p_to_keep)])
points_cached_series = np.empty([len(self.time_indices), len(p_to_keep),3])

for idx, tti in enumerate(self.time_indices): # oldest first
if idx > 0:
x_original_order, T_per_vertex = self.get_node_pos_and_temp(tti)
T_per_vertex_filt = [ T_per_vertex[i] for i in range(n_vertices) if i in p_to_keep ]
Temp_per_vertex_series[idx,:] = T_per_vertex_filt
point_original_to_cached = np.full(n_vertices,-1,dtype = np.int32)
count = 0
for i in range(n_vertices):
if (i in p_to_keep):
point_original_to_cached[i] = len(points_cached)
points_cached.append(x_original_order[i,:])
points_cached_series.append(np.array(points_cached))
points_cached_series[idx,count,:]=x_original_order[i,:]
point_original_to_cached[i]= count
count += 1


hexa_renumbered = [ [point_original_to_cached[i] for i in hexa] for hexa in hexa_to_keep ]

from os import path

filename_hex = path.join(out_path, self.modelName+'_hexa_ts_'+str(self.tti)+'.epc')
write_hexa_grid_with_timeseries(filename_hex, points_cached_series, hexa_renumbered, "hexamesh",
Temp_per_vertex_series )
Temp_per_vertex_series,age_per_vertex_keep, poro0_per_cell, decay_per_cell, density_per_cell,
cond_per_cell, rhp_per_cell, lid_per_cell )
return filename_hex

def heatflow_at_crust_sed_boundary(self):
Expand Down Expand Up @@ -1549,6 +1554,31 @@ def boundary(x):
return res.flatten()


def global_except_hook(exctype, value, traceback):
"""https://github.com/chainer/chainermn/issues/236
"""
try:
sys.stderr.write("\n*****************************************************\n")
sys.stderr.write("Uncaught exception was detected on rank {}. \n".format(
MPI.COMM_WORLD.Get_rank()))
from traceback import print_exception
print_exception(exctype, value, traceback)
sys.stderr.write("*****************************************************\n\n\n")
sys.stderr.write("\n")
sys.stderr.write("Calling MPI_Abort() to shut down MPI processes...\n")
sys.stderr.flush()
finally:
try:
MPI.COMM_WORLD.Abort(1)
except Exception as e:
sys.stderr.write("*****************************************************\n")
sys.stderr.write("Sorry, we failed to stop MPI, this process will hang.\n")
sys.stderr.write("*****************************************************\n")
sys.stderr.flush()
raise e

sys.excepthook = global_except_hook

def run_3d( builder:Builder, parameters:Parameters, start_time=182, end_time=0, pad_num_nodes=0,
out_dir = "out-mapA/",sedimentsOnly=False, writeout=True, base_flux=None):
logger.setLevel(10) # numeric level equals DEBUG
Expand Down Expand Up @@ -1612,9 +1642,9 @@ def run_3d( builder:Builder, parameters:Parameters, start_time=182, end_time=0,
comm.send(mm2.Tarr, dest=0, tag=((comm.rank-1)*10)+24)
if comm.rank==0:
mm2.receive_mpi_messages()
EPCfilename = mm2.write_hexa_mesh_resqml("temp/", end_time)
logger.info(f"RESQML model written to: {EPCfilename}")
# EPCfilename = mm2.write_hexa_mesh_resqml("temp/", end_time)
# logger.info(f"RESQML model written to: {EPCfilename}")
EPCfilename_ts = mm2.write_hexa_mesh_timeseries("temp/")
logger.info(f"RESQML partial model with timeseries written to: {EPCfilename_ts}")
read_mesh_resqml_hexa(EPCfilename) # test reading of the .epc file
read_mesh_resqml_hexa(EPCfilename_ts) # test reading of the .epc file
return mm2
108 changes: 92 additions & 16 deletions warmth/resqpy_helpers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
import resqpy.olio.uuid as bu
import resqpy.unstructured as rug
import resqpy.time_series as rts

from warmth.logging import logger
#
#
#
Expand Down Expand Up @@ -261,7 +261,8 @@ def read_mesh_resqml_hexa(epcfilename, meshTitle = 'hexamesh'):
# read properties
#

temp_uuid = model.uuid(title = 'Temperature')
temp_uuid = model.uuids(title = 'Temperature')
temp_uuid = temp_uuid[0]
assert temp_uuid is not None
temp_prop = rqp.Property(model, uuid = temp_uuid)
assert temp_prop.uom() == 'degC'
Expand All @@ -276,11 +277,12 @@ def read_mesh_resqml_hexa(epcfilename, meshTitle = 'hexamesh'):
assert layerID_prop.indexable_element() == 'cells'
print( layerID_prop.array_ref().shape, layerID_prop.array_ref()[0:10] ) # .array_ref() exposes the values as numpy array

titles=['Temperature', 'Age', 'LayerID', 'Porosity_initial', 'Porosity_decay', 'Density_solid', 'insulance_thermal', 'Radiogenic_heat_production']
for title in titles:
prop_uuid = model.uuid(title = title)
titles=[ 'Age', 'LayerID', 'Porosity_initial', 'Porosity_decay', 'Density_solid', 'insulance_thermal', 'Radiogenic_heat_production']
titles_uuid = [model.uuid(title = title) for title in titles]
titles_uuid.append(temp_uuid)
for prop_uuid in titles_uuid:
prop = rqp.Property(model, uuid = prop_uuid)
print(title, prop.indexable_element(), prop.uom(), prop.array_ref()[0:10] )
print(prop.title, prop.indexable_element(), prop.uom(), prop.array_ref()[0:10] )


def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexamesh",
Expand Down Expand Up @@ -459,7 +461,8 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame


def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle = "hexamesh",
Temp_per_vertex_series=None ):
Temp_per_vertex_series=None, age_per_vertex=None, poro0_per_cell=None, decay_per_cell=None, density_per_cell=None,
cond_per_cell=None, rhp_per_cell=None, lid_per_cell=None ):
"""Writes the given hexahedral mesh, defined by arrays of nodes and cell indices, into a RESQML .epc file
Given SubsHeat properties are optionally written.
Expand All @@ -477,8 +480,8 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
NOTE: writing properties that are defines per-node (have 'nodes' as indexable element) requires a patched version of resqpy!
"""

nodes = nodes_series[0]
node_count = len(nodes)
nodes_time_0 = nodes_series[-1,:,:] # present-day at last index
node_count = nodes_time_0.shape[0]
faces_per_cell = []
nodes_per_face = []
faces_dict = {}
Expand Down Expand Up @@ -514,8 +517,10 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
model = rq.new_model(filename)
crs = rqc.Crs(model)
crs.create_xml()
# present-day is set as 1 due resqpy
million_years_offset = 0
times_in_years = [ int(max((t+million_years_offset)*1e6,1+million_years_offset)) for t in list(range(nodes_series.shape[0]-1,-1,-1))]

times_in_years = [ max(int(t*1e6),1) for t in list(range(len(nodes_series)-1,-1,-1))]
gts = rts.GeologicTimeSeries.from_year_list(model, times_in_years, title="warmth simulation")
gts.create_xml()
rts.timeframe_for_time_series_uuid(model, gts.uuid)
Expand Down Expand Up @@ -543,7 +548,7 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
hexa.cell_face_is_right_handed = cell_face_is_right_handed # False for all faces for external cells

# points
hexa.points_cached = nodes
hexa.points_cached = nodes_time_0

# basic validity check
hexa.check_hexahedral()
Expand All @@ -556,12 +561,13 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
pc = hexa.property_collection

# nodes0 = nodes.copy()
for time_index in range(len(nodes_series)-1,-1,-1):
# nodes2 = nodes0 + [0,0,time_index*10]
nodes2 = nodes_series[time_index].astype(np.float32)
for time_index in range(nodes_series.shape[0]-1,-1,-1): #oldest first

nodes2 = nodes_series[time_index,:,:].astype(np.float32)

pc.add_cached_array_to_imported_list(nodes2,
'dynamic nodes',
'points',
"points",
uom = 'm',
property_kind = 'length',
realization = 0,
Expand All @@ -572,7 +578,7 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
tt = Temp_per_vertex_series[time_index].astype(np.float32)
pc.add_cached_array_to_imported_list(tt,
'Temperature',
'Temperature',
"Temperature",
uom = 'degC',
property_kind = 'thermodynamic temperature',
realization = 0,
Expand All @@ -582,5 +588,75 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
pc.write_hdf5_for_imported_list()
pc.create_xml_for_imported_list_and_add_parts_to_model(time_series_uuid = gts.uuid)


if age_per_vertex is not None:
_ = rqp.Property.from_array(model,
age_per_vertex.astype(np.float32),
source_info = 'SubsHeat',
keyword = 'Age',
support_uuid = hexa.uuid,
property_kind = 'geological age',
indexable_element = 'nodes',
uom = 'y')

if lid_per_cell is not None:
_ = rqp.Property.from_array(model,
lid_per_cell.astype(np.int32),
source_info = 'SubsHeat',
keyword = 'LayerID',
support_uuid = hexa.uuid,
property_kind = 'layer ID',
indexable_element = 'cells',
uom = 'Euc',
discrete=True)

if poro0_per_cell is not None:
_ = rqp.Property.from_array(model,
poro0_per_cell.astype(np.float32),
source_info = 'SubsHeat',
keyword = 'Porosity_initial',
support_uuid = hexa.uuid,
property_kind = 'porosity',
indexable_element = 'cells',
uom = 'm3/m3')
if decay_per_cell is not None:
_ = rqp.Property.from_array(model,
decay_per_cell.astype(np.float32),
source_info = 'SubsHeat',
keyword = 'Porosity_decay',
support_uuid = hexa.uuid,
property_kind = 'porosity decay',
indexable_element = 'cells',
uom = 'Euc')
if density_per_cell is not None:
_ = rqp.Property.from_array(model,
density_per_cell.astype(np.float32),
source_info = 'SubsHeat',
keyword = 'Density_solid',
support_uuid = hexa.uuid,
property_kind = 'density',
indexable_element = 'cells',
uom = 'kg/m3')
if cond_per_cell is not None:
#
# we write thermal conductivity as its inverse, the thermal insulance
#
_ = rqp.Property.from_array(model,
np.reciprocal(cond_per_cell).astype(np.float32),
source_info = 'SubsHeat',
keyword = 'insulance_thermal',
support_uuid = hexa.uuid,
property_kind = 'thermal insulance',
indexable_element = 'cells',
uom = 'deltaK.m2/W')
if rhp_per_cell is not None:
_ = rqp.Property.from_array(model,
rhp_per_cell.astype(np.float32),
source_info = 'SubsHeat',
keyword = 'Radiogenic_heat_production',
support_uuid = hexa.uuid,
property_kind = 'heat',
indexable_element = 'cells',
uom = 'W/m3')
model.store_epc()

Loading

0 comments on commit 93f8812

Please sign in to comment.