Skip to content
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

Vectorized mesh updates #17

Merged
merged 8 commits into from
Jan 17, 2024
Merged

Vectorized mesh updates #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
274f2f6
WIP: Fix use with MPI: mesh output working to EPC (no timeseries) wor…
Dec 27, 2023
60a9ada
Support re-uniting mesh and vertex properties from all MPI nodes in t…
Dec 27, 2023
2be8ec9
WIP: fix multi-timestep EPC output
Jan 3, 2024
5845ae5
Fix treatment of top and bottom vertex positions in MPI context.
Jan 4, 2024
d561d02
Fix output of timeseries EPC file in MPI multiprocess setting.
Jan 4, 2024
5b3a5e8
code cleanup
jtai-eq Jan 8, 2024
7cb41ad
support parallel execution in 1D simualor in MPI context; store 1D-mo…
jtai-eq Jan 16, 2024
cac5fbd
fix use of comm.barrier in the 1D simulations; pass use_mpi flag alr…
jtai-eq Jan 16, 2024
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
78 changes: 43 additions & 35 deletions tests/warmth3d/test_3d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,15 @@
import numpy as np
from mpi4py import MPI
import pickle
import time

def test_3d_compare():
comm = MPI.COMM_WORLD
if comm.rank == 0:
inc = 1000
model_pickled = f"model-out-inc_{inc}.p"
if comm.rank == 0 and not os.path.isfile(model_pickled):
global runtime_1D_sim
st = time.time()
maps_dir = Path("./docs/notebooks/data/")
model = warmth.Model()

Expand All @@ -22,8 +27,6 @@ def test_3d_compare():
inputs.loc[6]=[182,"182.gri",None,"Erosive"]
model.builder.input_horizons=inputs


inc = 1000
model.builder.define_geometry(maps_dir/"0.gri",xinc=inc,yinc=inc,fformat="irap_binary")

model.builder.extract_nodes(4,maps_dir)
Expand All @@ -50,12 +53,13 @@ def test_3d_compare():
model.parameters.tetha = 0
model.parameters.alphav = 0

model.simulator.run(save=True,purge=True, parallel=False)
model.simulator.run(save=True, purge=True, parallel=True, use_mpi=(comm.size>1))

import pickle
pickle.dump( model, open( "model-out.p", "wb" ) )
model = pickle.load( open( "model-out.p", "rb" ) )
runtime_1D_sim = time.time() - st
print("Total time 1D simulations:", runtime_1D_sim)

pickle.dump( model, open( model_pickled, "wb" ) )
model = pickle.load( open( model_pickled, "rb" ) )
try:
os.mkdir('mesh')
except FileExistsError:
Expand All @@ -66,41 +70,43 @@ def test_3d_compare():
pass

comm.Barrier()
import pickle
model = pickle.load( open( "model-out.p", "rb" ) )
model = pickle.load( open( model_pickled, "rb" ) )
comm.Barrier()
# if comm.rank == 0:
# mm2 =
mm2, posarr, Tarr = run_3d(model.builder,model.parameters,start_time=model.parameters.time_start,end_time=0, pad_num_nodes=2,writeout=False, base_flux=None)

mm2 = run_3d(model.builder,model.parameters,start_time=model.parameters.time_start,end_time=0, pad_num_nodes=2,writeout=False, base_flux=None)

nnx = (model.builder.grid.num_nodes_x+2*mm2.padX)
nny = (model.builder.grid.num_nodes_y+2*mm2.padX)
# hx = model.builder.grid.num_nodes_x // 2
# hy = model.builder.grid.num_nodes_y // 2
hx = nnx // 2
hy = nny // 2
comm.Barrier()
if comm.rank == 0:
nnx = (model.builder.grid.num_nodes_x+2*mm2.padX)
nny = (model.builder.grid.num_nodes_y+2*mm2.padX)
hx = nnx // 2
hy = nny // 2

nn = model.builder.nodes[hy-mm2.padX][hx-mm2.padX]
dd = nn._depth_out[:,0]

mm2_pos, mm2_temp = mm2.get_node_pos_and_temp(-1)

nn = model.builder.nodes[hy-mm2.padX][hx-mm2.padX]
dd = nn._depth_out[:,0]
node_ind = hy*nnx + hx
v_per_n = int(mm2_pos.shape[0]/(mm2.num_nodes))

node_ind = hy*nnx + hx
# v_per_n = int(mm2.mesh_vertices.shape[0]/(model.builder.grid.num_nodes_y*model.builder.grid.num_nodes_x))
v_per_n = int(mm2.mesh_vertices.shape[0]/(mm2.num_nodes))
temp_1d = np.nan_to_num(nn.temperature_out[:,0], nan=5.0)

# temp_3d_ind = np.array([ np.where([mm2.mesh_reindex==i])[1][0] for i in range(node_ind*v_per_n, (node_ind+1)*v_per_n) ] )
temp_3d_ind = np.array([ i for i in range(node_ind*v_per_n, (node_ind+1)*v_per_n) ] )

temp_1d = np.nan_to_num(nn.temperature_out[:,0], nan=5.0)
temp_3d_ind = np.array([ np.where([mm2.mesh_reindex==i])[1][0] for i in range(node_ind*v_per_n, (node_ind+1)*v_per_n) ] )
dd_mesh = mm2.mesh.geometry.x[temp_3d_ind,2] - mm2.sed_diff_z[range(node_ind*v_per_n, (node_ind+1)*v_per_n)]
temp_3d_mesh = mm2.u_n.x.array[temp_3d_ind]
# dd_mesh = mm2.mesh.geometry.x[temp_3d_ind,2] - mm2.sed_diff_z[range(node_ind*v_per_n, (node_ind+1)*v_per_n)]
dd_mesh = mm2_pos[temp_3d_ind,2] - mm2.sed_diff_z[range(node_ind*v_per_n, (node_ind+1)*v_per_n)]
# temp_3d_mesh = mm2.u_n.x.array[temp_3d_ind]
temp_3d_mesh = mm2_temp[temp_3d_ind]

temp_1d_at_mesh_pos = np.interp(dd_mesh, dd, temp_1d)
dd_subset = np.where(dd_mesh<5000)
temp_1d_at_mesh_pos = np.interp(dd_mesh, dd, temp_1d)
dd_subset = np.where(dd_mesh<5000)

max_abs_error = np.amax(np.abs(temp_1d_at_mesh_pos-temp_3d_mesh))
max_abs_error_shallow = np.amax(np.abs(temp_1d_at_mesh_pos[dd_subset]-temp_3d_mesh[dd_subset]))
max_abs_error = np.amax(np.abs(temp_1d_at_mesh_pos-temp_3d_mesh))
max_abs_error_shallow = np.amax(np.abs(temp_1d_at_mesh_pos[dd_subset]-temp_3d_mesh[dd_subset]))

print("Max error overall:", max_abs_error)
print("Max error <5000m:", max_abs_error_shallow)
print("Max error overall:", max_abs_error)
print("Max error <5000m:", max_abs_error_shallow)

# assert (max_abs_error<25.0), "Temperature difference between 3D and 1D simulations is >25"
# assert (max_abs_error_shallow<5.0), "Temperature difference between 3D and 1D simulations is >5 in the sediments."
Expand All @@ -109,4 +115,6 @@ def test_3d_compare():


if __name__ == "__main__":
test_3d_compare()
test_3d_compare()
if 'runtime_1D_sim' in globals():
print("Total time 1D simulations:", runtime_1D_sim)
Loading
Loading