import time
import logging
from pathlib import Path
from typing import Union
from math import log10, ceil
import numpy as np
import pygame
# from wolfhece.PyGui import Wolf2DModel # deprecated in favor of prev_sim2D
from wolfhece.mesh2d.wolf2dprev import prev_sim2D
from importlib.metadata import version
[docs]
wolf_hece_version = tuple(map(int, version("wolfhece").split('.')))
if wolf_hece_version <= (1,8,11):
from wolfhece.mesh2d.wolf2dprev import prev_boundary_condition
elif wolf_hece_version >= (2,0,5):
from wolfhece.mesh2d.wolf2dprev import boundary_condition_2D as prev_boundary_condition
else:
# 1.8.13, 2.0.1 are broken for us :-(
raise Exception(f"I don't recognize this version of wolfhece : {wolf_hece_version}")
from ..glsimulation import GLSimulation, get_report_frequency
from ..tile_packer import TilePackingMode
from ..loaders.utils import convert_bc_type, make_impervious_bathymetry
from ..simple_simulation import InfiltrationInterpolation, SimulationDuration, SimulationDurationType, TimeStepStrategy, BoundaryConditionsTypes
[docs]
def load_sim_to_gpu(
model: Union[Path, prev_sim2D],
init_gl=True,
glsim: Union[None, GLSimulation] = None,
tile_size=16,
interp_infiltration: InfiltrationInterpolation = InfiltrationInterpolation.NONE,
shader_log_path=None,
interpret_freq_type_as_wolf: bool = False,
tiles_packing_mode: TilePackingMode = TilePackingMode.REGULAR,
optimize_indirection: bool = False,
fail_if_invalid_sim: bool = True
) -> GLSimulation:
"""
`model`: You can pass either a directory to a model (a `Path`) or
a `Wolf2DModel`.
`glsim`: If provided, the model will be loaded into the given
simulation. If not provided, a brand new simulation will be
created.
`interpret_freq_type_as_wolf`: If true and unit of report frequency is
seconds (instead of steps), then we interpret the duration of the simulation
as a number of seconds (instead of a number of steps). This param was introduced
because in our test/debug, "npas" is always expected to be a number of steps
(and not seconds).
"""
assert isinstance(model, Path) or isinstance(model, prev_sim2D)
assert interp_infiltration in (InfiltrationInterpolation.NONE, InfiltrationInterpolation.LINEAR)
errors = []
def check_error(cond, msg):
if not cond:
errors.append(msg)
if init_gl:
pygame.init()
# Loads a Wolf simulation into the simulator
start_time = time.time()
if isinstance(model, Path):
logging.info(f"Loading model data from Wolf model {model}")
# Will figure out itself the name of the files inside the directory
wolf2d=prev_sim2D(fname=str(model))
#print(wolf2d.parameters.to_yaml())
# This is hackish but allows to load
# the data of the simultaion (seems
# like things are lazy loaded, or something like
# that)
else:
wolf2d = model
wolf2d.force_load()
logging.info(f"Uploading to GPU")
# TODO At some point, I'll have to merge all
# of that back into PA's models, this is s huge duplication
# of code (it was done to better understand what I take out
# of PA's model; now that I understand them better, I could
# directly take data out of them)
def copy_and_extend(a):
plus_row = np.vstack( [ a, np.zeros( (1,a.shape[1]) ) ])
plus_col = np.hstack( [ plus_row, np.zeros( (plus_row.shape[0],1) ) ])
return plus_col
if wolf2d.nb_blocks > 1:
logging.warning("The modelisation appears to be multi-blocks. I will *ONLY* handle the *first* one !")
# FIXME Depending on the source, you may or may not
# have a valid maks here.
nap = wolf2d.napbin.array.data.copy()
# FIXME Depending on the source, an active (= non masked) value is denoted
# by +1 or -1; why are there two possibilities ? Inactive
# seem to be always zero though.
# So we fix that: 0 means inactive; 1 means active.
nap[nap != 0] = 1
check_error(not np.all(nap == 0), "NAP matrix is full of zeros, meaning every cells are inactive!")
# FIXME Why copy ?
bathymetry = wolf2d.top.array.filled( np.max(wolf2d.top.array) )
check_error(np.sum(bathymetry < 0) == 0, "I have found negative levels of bathymetry. I can't handle that.")
water_height = wolf2d.hbin.array.filled(0)
# print(water_height)
# print(wolf2d.hbin.array)
check_error( np.sum(water_height < 0) == 0, "I have found negative levels of water. I can't handle that.")
qxbin = wolf2d.qxbin.array.filled(0)
qybin = wolf2d.qybin.array.filled(0)
manning = wolf2d.frot.array.filled( np.max(wolf2d.frot.array) )
check_error( np.sum(manning < 0) == 0, "I have found negative Manning coefficients. I can't handle that.")
# wolf2d.infiltration.plot_plt()
infiltration = wolf2d.infiltration
noptpas=wolf2d.parameters._scheme_optimize_timestep
nxfin=wolf2d.nbx
nyfin=wolf2d.nby
dxfin=float(wolf2d.dx)
dyfin=float(wolf2d.dy)
check_error( dxfin == dyfin, "Only squared meshes supported; 'dxfin' must equal 'dyfin'")
ponderation = wolf2d.parameters._scheme_rk
froude_max = wolf2d[1]._froude_max
if ponderation > 1.0 or ponderation <= 0.0:
logging.warning(f"Unsupported ponderation value {ponderation}. I expect 'ponderation' between zero and one. I'll replace by R-K with 0.3/0.7.")
ponderation = 0.3
vncsouhaite=wolf2d.parameters._scheme_cfl
npas= int(wolf2d.parameters._nb_timesteps)
tstep = wolf2d.parameters._timestep_duration or None
# This mumbo jumbo to make sure wxWidget leaves us in peace to define our
# own GL context It doesn't work if we can't close all frames which happens
# when your mouse is on an OpenGL frame and the text lens gets activated. So
# keep your mouse in the corner of the screen (or find a way to close all
# frames)
# wolf2d.allviews.Close()
# wolf2d.allviews.mytooltip.Close()
# wolf2d.Close()
# wolf2d.mylogs.GetFrame().Close()
# ex = None
# wolf2d = None
# logging.info("Sleeping")
# # sleep(5)
# logging.info("Awaking")
if init_gl:
from math import sqrt
display_nfo = pygame.display.Info()
screen_width = int(display_nfo.current_w*0.9)
screen_height = int(display_nfo.current_h*0.9)
sim_area = nxfin*nyfin
screen_area = (screen_width * screen_height)/5
area_ratio = screen_area / sim_area
aspect_ratio = nxfin / nyfin
window_height = sqrt(area_ratio)*nyfin
window_width = sqrt(area_ratio)*aspect_ratio*nyfin
if window_height > screen_height:
window_height = screen_height
window_width = aspect_ratio * window_height
elif window_width > screen_width:
window_width = screen_width
window_height= window_width / aspect_ratio
window_height = int(round(window_height))
window_width = int(round(window_width))
# w,h=nxfin,nyfin
# |pygame.DOUBLEBUF
# Note: vsync=0: not sure this works. But setting no vblank on the
# nvidia control panel does work.
#window_width, window_height = nxfin, nyfin
pygame.display.set_mode((window_width, window_height), pygame.OPENGL|pygame.DOUBLEBUF, vsync = 0)
#pygame.display.set_mode((w, h), pygame.OPENGL)
# # The glDisables will work only after pygame.display.set_mode has been
# # called.
# glDisable(GL_BLEND)
# glDisable(GL_DEPTH_TEST)
# # We absolutely don't want anti alisaing (for exampel
# # when computeing the max step size)
# glDisable(GL_POLYGON_SMOOTH)
# logging.info("OpenGL initialized")
# numpy is row-major (https://ncar-hackathons.github.io/scientific-computing/numpy/02_memory_layout.html)
# consequently, addressing a numpy array is done as a[y,x] (y first).
# In a WolfArray, the x and y are swapped when compared
# to the width and height of a numpy array.
# I have two solutions:
# - I keep the (x,y) shape but then I need to transpose all WolfArray matrices
# before sending the to GPU wolf and transpose them again when writing
# results.
# - I swap x,y to (y,x) but then I am at the mercy of any non-symetric
# process in the simulation.
check_error( (nxfin, nyfin) == bathymetry.shape, f"Shapes differ ? Params say: {wolf2d.nby}x{wolf2d.nbx}, .top file says {bathymetry.shape[0]}x{bathymetry.shape[1]}")
# FIXME contiguous array ? Not sure it is of anyhelp since
# pyopengl inputs np.array.
nap = np.ascontiguousarray(nap.transpose())
logging.info(f"NAP as numpy array is {nap.shape[0]} rows x {nap.shape[1]} columns")
if False:
# debug code
import matplotlib.pyplot as plt
print(f"{np.min(nap.data)} - {np.max(nap.data)}")
plt.title("Showing NAP values to figure which is which")
plt.imshow(nap, origin="lower")
plt.colorbar()
plt.show()
bathymetry = np.ascontiguousarray(bathymetry.transpose())
if False:
numpy_crd = (4093-1,631-1)
b = bathymetry[numpy_crd]
from matplotlib import pyplot as plt
plt.imshow(np.flipud(bathymetry))
plt.show()
water_height = np.ascontiguousarray(water_height.transpose())
qxbin = np.ascontiguousarray(qxbin.transpose())
qybin = np.ascontiguousarray(qybin.transpose())
manning = np.ascontiguousarray(manning.transpose())
# Replace "null" value (-99999) byt something
# meaningful so that when we query cells outside
# the "general contour", we find somethign appropriate.
# FIXME This should be done in the shaders so that computations
# are adapted to missing data (instead of providing "default"
# data). Note we have discussed that with PA in dec/2023 and
# we decided it was OK for the moment.
# Create a wall that is ten times has tall as the max bathy.
# Additional sorcery to make that height something clearly
# artificial, ie 9999.
#unreachable_height = (10 ** (ceil(log10(np.max( bathymetry))) + 1)) - 1
#unreachable_height = round(np.max( bathymetry) + 2*np.max(water_height))
# bathymetry[bathymetry < 0] = unreachable_height
# water_height[water_height < 0] = 0
# qxbin[qxbin == -99_999.0] = 0
# qybin[qybin == -99_999.0] = 0
# When handling boundary condition, we may look in cells
# outside the computation domain (inactive, masked cells).
# For these cells, we expect
# "neutral" data (zeros) that won't affect computations.
msk = nap == 0 # => True = inactive cell
water_height[msk] = 0
qxbin[msk] = 0 # BUG ???
qybin[msk] = 0
# We enclose the simulation domain in a border
# Make sure to leave the B.C. as they are.
# FIXME I don't think it is necessary...
nap_draft = np.copy(nap)
for bc in wolf2d.parameters.clfbx.mybc:
cell_numpy_crd = (bc.j-1, bc.i-1)
nap_draft[cell_numpy_crd] = 1
for bc in wolf2d.parameters.clfby.mybc:
cell_numpy_crd = (bc.j-1, bc.i-1)
nap_draft[cell_numpy_crd] = 1
bathymetry = make_impervious_bathymetry(nap_draft, bathymetry)
# min_pos = (628-1, 1887-1)
# sl_row2 = slice(min_pos[0]-4, min_pos[0]+5)
# sl_col2 = slice(min_pos[1]-4, min_pos[1]+5)
# dbg_bath = bathymetry[sl_row2, sl_col2]
# dbg_qxbin = qxbin[sl_row2, sl_col2]
# dbg_qybin = qybin[sl_row2, sl_col2]
# dbg_h = water_height[sl_row2, sl_col2]
# These two are important because Qx is taken into account
# when doing splitting. So if we need to split between
# a cell that is inside the domain and one that is outside,
# having zeros in the (outside domain) Qx,Qy makes computation
# easier
if not np.all( np.logical_not( (qxbin != 0) & (water_height == 0) )):
logging.warning("There are cells where Qx is non-zero but where's no water ! Here are some coordinates:")
logging.warning(np.fliplr(np.argwhere((qxbin != 0) & (water_height == 0))))
qxbin[(qxbin != 0) & (water_height == 0)] = 0
if not np.all( np.logical_not( (qybin != 0) & (water_height == 0) )):
logging.warning("There are cells where Qy is non-zero but where's no water ! Here are some coordinates:")
logging.warning(np.fliplr(np.argwhere((qybin != 0) & (water_height == 0))))
qybin[(qybin != 0) & (water_height == 0)] = 0
logging.info(f"{100*(np.sum(nap > 0)/(nxfin*nyfin)):.1f}% active meshes; {100*(np.sum(water_height !=0)/(nxfin*nyfin)):.1f}% meshes with water")
# assert np.alltrue( np.logical_not( (qxbin != 0) & (water_height == 0) )), \
# "There are cells where Qx is non-zero but where's no water !"
# assert np.alltrue( np.logical_not( (qybin != 0) & (water_height == 0) )), \
# "There are cells where Qy is non-zero but where's no water !"
# Handle weak boundary conditions, horizontal direction
# 1=active cell, shall be computed; 0=inactive cell, will not be computed.
# Pay attention to the fact that if an active cell touches an inactive
# one, then, there must be boundary conditions.
# Bit 1: cell is active(==1) or not (==0)
# Bit 2: weak boundary condition QX (left border) (==1) or not (==0)
weak_boundary_cond = np.zeros((nyfin, nxfin), dtype=np.uint32)
weak_boundary_cond[(water_height > 0) | (nap == 1)] = 1 # FIXME Not correct, use the NAP mask !
# Mapping cell numpy coordinates (x,y) onto an array
# of values. Each value represents a boundary condition.
# The index of the value represents the type of the condition.
# If the value is BC_BLANK_VALUE it means the b.c. is not
# used on the cell.
# e.g. bc_descriptions[ (i,j) ] = [BC's H on left border's value, HMOD on left border's value, ... ]
bc_descriptions = dict()
BC_BLANK_VALUE=np.float32(-9_999_999) # Exactly represented byt float32
NB_BC_TYPES=10
BC_TABLE_INDEX_FOR_H_ON_LEFT = 1
BC_TABLE_INDEX_FOR_HMOD_ON_LEFT = 2
BC_TABLE_INDEX_FOR_QX_ON_LEFT = 0
BC_TABLE_INDEX_FOR_QY_ON_LEFT = 3 # Qy on a Y horizontal border FIXME we're missing it on the vertical border !
BC_TABLE_INDEX_FOR_FROUDE_NORMAL_ON_LEFT = 8
BC_TABLE_INDEX_FOR_H_ON_BOTTOM = 4
BC_TABLE_INDEX_FOR_HMOD_ON_BOTTOM = 5
BC_TABLE_INDEX_FOR_QX_ON_BOTTOM = 6 # Qx applied to a bottom border (so Ax parallel to the border)
BC_TABLE_INDEX_FOR_QY_ON_BOTTOM = 7 # Qy applied to a bottom border (so Ax parallel to the border)
BC_TABLE_INDEX_FOR_FROUDE_NORMAL_ON_BOTTOM = 9
# ------------------------------------------------------------------------
# Boundary Conditions on left borders
# ------------------------------------------------------------------------
bc: prev_boundary_condition
for bc in wolf2d.parameters.clfbx.mybc:
bc_type = convert_bc_type(bc.ntype)
logging.debug(f"B.C. left border at {bc.i},{bc.j}, type: {bc_type}")
# Attention! x,y are swapped !!!
cell_numpy_crd = (bc.j-1, bc.i-1)
if cell_numpy_crd not in bc_descriptions:
bc_descriptions[cell_numpy_crd] = [BC_BLANK_VALUE] * NB_BC_TYPES
x, y = bc.i-1, bc.j-1 # i,j are Fortran coords => translate to numpy
# - Either the BC is on the left border of a cell which
# is on the interior side of the border of the domain.
# - Either the BC is on the left border of a cell which
# is on the exterior side of the border of the domain.
left_cell_crd = (y, x - 1)
if not( (x == 0 and nap[cell_numpy_crd] == 1) \
or (nap[cell_numpy_crd] == 1 and nap[left_cell_crd] == 0) \
or (nap[cell_numpy_crd] == 0 and nap[left_cell_crd] == 1)):
err = f"Weak boundary condition at i={bc.i},j={bc.j}: this B.C. is not on the border of the computation domain"
logging.error(err)
raise Exception(err)
if bc_type == BoundaryConditionsTypes.QX:
logging.debug(f"Weak boundary condition over QX at i={bc.i},j={bc.j}, value={bc.val}")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_QX_ON_LEFT] = bc.val
elif bc_type == BoundaryConditionsTypes.QY:
logging.debug(f"Weak boundary condition over QY at i={bc.i},j={bc.j}, value={bc.val}")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_QY_ON_LEFT] = bc.val
elif bc_type == BoundaryConditionsTypes.H:
loc_msg = f"Weak boundary condition over H at i={bc.i},j={bc.j}"
logging.debug(f"{loc_msg}, value={bc.val}")
if bc.val < 0:
# negative values are a hack to express impervious boundaries
logging.error(f"{loc_msg}, h={bc.val} is negative !!!")
if qxbin[cell_numpy_crd] > 0:
logging.warning(f"{loc_msg} doesn't provide enough information to the cell (see cell flow direction)")
if nap[cell_numpy_crd] == 0 \
and bathymetry[left_cell_crd] > bc.val:
logging.error(f"{loc_msg}: seems located rightside of the " \
"domain but the bathymetry in the neighbour " \
f"cell is above the surface level of the b.c. {bathymetry[left_cell_crd]} > {bc.val}, I expect it to be below.")
if nap[cell_numpy_crd] == 1 \
and bathymetry[cell_numpy_crd] > bc.val:
logging.error(f"{loc_msg}: seems located leftside of the " \
"domain but the bathymetry in the cell is " \
f"above the surface level of the b.c. {bathymetry[left_cell_crd]} > {bc.val}, I expect it to be below.")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_H_ON_LEFT] = bc.val
elif bc_type == BoundaryConditionsTypes.HMOD:
logging.debug(f"Weak boundary condition over HMOD at i={bc.i},j={bc.j}, value={bc.val}")
if bc.val < 0:
# negative values are a hack to express impervious boundaries
logging.warning(f"Weak boundary condition over HMOD at i={bc.i},j={bc.j}, value={bc.val} is negative !!!")
if x < qxbin.shape[1] and y < qxbin.shape[0] and qxbin[y,x] > 0:
logging.warning(f"Weak boundary condition over HMOD at i={bc.i},j={bc.j} doesn't provide enough information to the cell (see cell flow direction)")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_HMOD_ON_LEFT] = bc.val
elif bc_type == BoundaryConditionsTypes.FROUDE_NORMAL:
logging.debug(f"Weak boundary condition Froude normal direction at i={bc.i},j={bc.j}, value={bc.val}")
check_error( bc.val > 0, f"The Froude number must be > 0, you gave {bc.val}")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_FROUDE_NORMAL_ON_LEFT] = bc.val
else:
logging.error(f"Unsupported weak condition type {bc_type}")
# Attention: if we have a BC on the left border of cell (i,j)
# we record it. But we also mark the cell (i-1,j) with a
# no-values boundary condition. This special "no-values b.c." cell
# will be interpreted by the GPU engine as a cell with a b.c.
# on its *right* border. This is a bit hacky but makes sure
# that the values of a boundary condition are not duplicated
# across cells (which always carries a risk of mis-duplication bug)
# Of course if one wants to change the BC on the (i-1,j) cell
# he can as it will not disable the detection on a right border BC.
# FIXME This is dangerous as usually a BC is defined on a border
# but is meant to produce effect only on one side of the border.
cell_numpy_crd_left = (cell_numpy_crd[0], cell_numpy_crd[1]-1)
if cell_numpy_crd_left not in bc_descriptions:
bc_descriptions[cell_numpy_crd_left] = [BC_BLANK_VALUE] * NB_BC_TYPES
# ------------------------------------------------------------------------
# Boundary Conditions on bottom borders
# ------------------------------------------------------------------------
bc: prev_boundary_condition
for bc in wolf2d.parameters.clfby.mybc:
bc_type = convert_bc_type(bc.ntype)
logging.debug(f"B.C. bottom border at {bc.i},{bc.j}, type: {bc_type}")
# Attention! x,y are swapped !!!
cell_numpy_crd = (bc.j-1, bc.i-1)
if cell_numpy_crd not in bc_descriptions:
bc_descriptions[cell_numpy_crd] = [BC_BLANK_VALUE] * NB_BC_TYPES
x, y = bc.i-1, bc.j-1 # i,j are Fortran coords => translate to numpy
bottom_cell_crd = (y - 1, x)
# - Either the BC is on the bottom border of a cell which
# is on the interior side of the border of the domain.
# - Either the BC is on the bottom border of a cell which
# is on the exterior side of the border of the domain.
if not( (x == 0 and nap[cell_numpy_crd] == 1) \
or (nap[cell_numpy_crd] == 1 and nap[bottom_cell_crd] == 0) \
or (nap[cell_numpy_crd] == 0 and nap[bottom_cell_crd] == 1)):
err = f"Weak boundary condition at i={bc.i},j={bc.j}: this B.C. is not on the border of the computation domain"
logging.error(err)
raise Exception(err)
if bc_type == BoundaryConditionsTypes.QX:
logging.debug(f"Weak boundary condition over QX (bottom border) at i,j:{(bc.i, bc.j)}, value={bc.val}")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_QX_ON_BOTTOM] = bc.val
elif bc_type == BoundaryConditionsTypes.QY:
logging.debug(f"Weak boundary condition over QY (bottom border) at i,j:{(bc.i, bc.j)}, value={bc.val}")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_QY_ON_BOTTOM] = bc.val
# elif bc_type == BoundaryConditionsTypes.H.value:
# logging.info(f"Weak boundary condition over H (bottom border) at i,j:{(bc.i, bc.j)}, value={bc.val}")
# if bc.val < 0:
# # negative values are a hack to express impervious boundaries
# logging.warning(f"Weak boundary condition over H (horiz border) at i={bc.i},j={bc.j}, value={bc.val} is negative !!!")
# if qxbin[y,x] > 0:
# logging.warning(f"Weak boundary condition over H (horiz border) at i={bc.i},j={bc.j} doesn't provide enough information to the cell (see cell flow direction)")
# if bathymetry[cell_numpy_crd] > bc.val:
# logging.error(f"Weak boundary condition over H at i={bc.i},j={bc.j}: the bathymetry is above the surface level, I expect it to be below.")
# bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_H_ON_BOTTOM] = bc.val
elif bc_type == BoundaryConditionsTypes.H:
# ........................................................................ NEW
loc_msg = f"Weak horizontal boundary condition over H (surface altitude) at i={bc.i},j={bc.j}"
#logging.info(f"{loc_msg}, value={bc.val}")
if bc.val < 0:
# negative values are a hack to express impervious boundaries
logging.error(f"{loc_msg}, h={bc.val} is negative !!!")
if qybin[cell_numpy_crd] > 0:
logging.warning(f"{loc_msg} doesn't provide enough information to the cell (see cell flow direction)")
if nap[cell_numpy_crd] == 0 \
and bathymetry[bottom_cell_crd] > bc.val:
logging.error(f"{loc_msg}: seems located top side of a mesh in the " \
"domain but the bathymetry in the neighbour " \
f"cell is above the surface level of the b.c. {bathymetry[bottom_cell_crd]} > {bc.val}, I expect it to be below.")
if nap[cell_numpy_crd] == 1 \
and bathymetry[cell_numpy_crd] > bc.val:
logging.error(f"{loc_msg}: seems located bottom side of a mesh in the " \
"domain but the bathymetry in the cell is " \
f"above the surface level of the b.c. {bathymetry[bottom_cell_crd]} > {bc.val}, I expect it to be below.")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_H_ON_BOTTOM] = bc.val
# ........................................................................ NEW
elif bc_type == BoundaryConditionsTypes.HMOD:
logging.debug(f"Weak boundary condition over HMOD (water height, horiz border) at i={bc.i},j={bc.j}, value={bc.val}")
if bc.val < 0:
# negative values are a hack to express impervious boundaries
logging.warning(f"But at i={bc.i},j={bc.j}, B.C.'s value={bc.val} is negative !!!")
if x < qybin.shape[1] and y < qybin.shape[0] and qybin[y,x] > 0:
logging.warning(f"But at i={bc.i},j={bc.j} B.C. doesn't provide enough information to the cell (see cell flow direction)")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_HMOD_ON_BOTTOM] = bc.val
elif bc_type == BoundaryConditionsTypes.FROUDE_NORMAL:
logging.debug(f"Weak boundary condition Froude normal direction at i={bc.i},j={bc.j}, value={bc.val}")
check_error( bc.val > 0, f"The Froude number must be > 0, you gave {bc.val}")
bc_descriptions[cell_numpy_crd][BC_TABLE_INDEX_FOR_FROUDE_NORMAL_ON_BOTTOM] = bc.val
else:
raise Exception(f"Unsupported weak condition type {bc_type} at i={bc.i},j={bc.j}")
# Mirror this BC on the top border of the bottom side cell.
cell_numpy_crd_bottom = (cell_numpy_crd[0]-1, cell_numpy_crd[1])
if cell_numpy_crd_bottom not in bc_descriptions:
bc_descriptions[cell_numpy_crd_bottom] = [BC_BLANK_VALUE] * NB_BC_TYPES
# ........................................................................
# Handle strong boundary conditions
# ........................................................................
strong_boundary_cond = np.ones((nyfin, nxfin), dtype=np.float32)
strong_boundary_cond[:,:] = -99999.0
strong_boundary_cond_values = np.zeros((nyfin, nxfin, 3), dtype=np.float32)
bc: prev_boundary_condition
for bc in wolf2d.parameters.clf.mybc:
bc_type = convert_bc_type(bc.ntype)
if bc_type == BoundaryConditionsTypes.QX:
x, y = bc.i-1, bc.j-1
logging.error(f"IGNORED Strong boundary condition over QX at x={x},y={y}={bc.val}")
strong_boundary_cond[x, y] = 0.0
strong_boundary_cond_values[y, x] = np.array([0.0,bc.val,0.0])
elif bc_type == BoundaryConditionsTypes.H:
# FIXME BoundaryConditionsTypes.H.value is shared with *weak*
# boundary conditions. So is that number OK to be shared ?
x, y = bc.i-1, bc.j-1
logging.error(f"IGNORED Strong boundary condition over H at x={x},y={y}={bc.val}")
strong_boundary_cond[x, y] = 1.0
strong_boundary_cond_values[y, x] = np.array([bc.val,0.0,0.0])
else:
raise Exception(f"Unsupported strong condition type {bc_type}")
# ........................................................................
# Checking boundary conditions are coherent
# ........................................................................
for crd, bc_values in bc_descriptions.items():
if np.all( np.array((bc_values) == BC_BLANK_VALUE)):
# We don't have any b.c. here.
# Two possible reasons:
# - Wolf has more types of BC than what we support here). So none
# is marked
# - this `bc_values` represents a beacon for a top or right b.c.
continue
ij_crd = "i:{},j:{}".format(crd[1]+1,crd[0]+1)
# Either there's no Qx and no Qy, either there is a Qx AND a Qy
check_error( (bc_values[BC_TABLE_INDEX_FOR_QX_ON_LEFT] != BC_BLANK_VALUE) == (bc_values[BC_TABLE_INDEX_FOR_QY_ON_LEFT] != BC_BLANK_VALUE), \
f"At {ij_crd}, If you set Qx on the left border you must set Qy too, and vice versa.")
check_error( (bc_values[BC_TABLE_INDEX_FOR_QX_ON_BOTTOM] != BC_BLANK_VALUE) == (bc_values[BC_TABLE_INDEX_FOR_QY_ON_BOTTOM] != BC_BLANK_VALUE), \
f"At {ij_crd}, If you set Qx (or Qy) on the bottom border you must set Qy (or Qx) too.")
h_or_hmod = bc_values[BC_TABLE_INDEX_FOR_H_ON_LEFT] != BC_BLANK_VALUE or bc_values[BC_TABLE_INDEX_FOR_HMOD_ON_LEFT] != BC_BLANK_VALUE
qx_or_qy = bc_values[BC_TABLE_INDEX_FOR_QX_ON_LEFT] != BC_BLANK_VALUE or bc_values[BC_TABLE_INDEX_FOR_QY_ON_LEFT] != BC_BLANK_VALUE
# Read : (none of them) OR (any of them but not the two to of them)
check_error( (not h_or_hmod and not qx_or_qy) or (h_or_hmod ^ qx_or_qy), \
f"At {ij_crd} {bc_values}, If you set Qx,Qy on left border then you can't set h/hmod.")
if bc_values[BC_TABLE_INDEX_FOR_QX_ON_LEFT] != BC_BLANK_VALUE:
if nap[crd] != 0:
check_error( nap[crd[0],crd[1]-1] == 0, \
f"Setting a left border b.c. at {ij_crd} makes sense if leftside cell is out of computation domain")
else:
check_error( nap[crd[0],crd[1]-1] != 0, \
f"Setting a right border b.c. at {ij_crd} makes sense if leftside cell is inside of computation domain")
# ........................................................................
# Converting boundary conditions into GPU data
# ........................................................................
# Convert the noted cells' boundary conditions into a
# sequential array and pointing towards it.
if bc_descriptions:
bc_cells_array = []
for i,(cell_numpy_crd,bc_values) in enumerate(bc_descriptions.items()):
bc_cells_array.append(bc_values)
# Here is a tricky thing:
# If there are some b.c. on a cell, then we give it an index >= 1.
# If there are no b.c. we give the index 0.
# This way we can know if there are b.c. or not without
# looking at b.c. themselves. It's an optimisation.
# Moreover, don't forget that a right-border B.C. will
# be located outside of the active cells !
weak_boundary_cond[cell_numpy_crd] |= (i+1) << 1 # bit 1 is NAP.
logging.debug(f"weak b.c.: #{weak_boundary_cond[cell_numpy_crd]}: {cell_numpy_crd} -> {bc_values}")
bc_cells_array = np.array(bc_cells_array,dtype=np.float32)
#print(bc_cells_array)
#bc_cells_array = bc_cells_array.transpose()
logging.info(f"There are {len(bc_descriptions)} boundary conditions.")
else:
bc_cells_array = None
logging.debug(f"{100*(np.sum(weak_boundary_cond == 1)/(nxfin*nyfin)):.1f}% active meshes without BC in weak boundary conditions array")
if errors:
for error in errors:
logging.error(error)
if errors and fail_if_invalid_sim:
raise Exception("Errors detected, check log")
# ........................................................................
# Infiltration
# ........................................................................
if hasattr(infiltration, "my_Q") and infiltration.my_Q.size > 0:
infiltration_timings = infiltration.my_Q
nb_infiltration_lines = infiltration_timings.shape[0]
logging.debug("Infiltration timings array")
logging.debug(infiltration_timings)
assert np.all(np.diff(infiltration_timings[:,0]) > 0), f"The infiltration timings are not strictly increasing: {infiltration_timings[:,0]}"
if interp_infiltration == InfiltrationInterpolation.LINEAR:
assert infiltration_timings[0,0] == 0, f"The first infiltration time must be zero if you use interpolation"
assert nb_infiltration_lines >= 2, f"You need at least two data points in infiltrations if you use interpolation"
zones_numbers = np.unique(wolf2d.inf.array[wolf2d.inf.array > 0]).tolist()
logging.info(f"Infitration zones numbers are : {zones_numbers}")
used_zones = set()
for row in infiltration_timings:
for zone_ndx in range(1, row.shape[0]):
if row[zone_ndx] > 0:
used_zones.add(zone_ndx)
if used_zones != set(zones_numbers):
logging.warning(f"Mismatch between the zones in the infiltration table ({','.join(map(str,used_zones))}) and the zones in the infiltration map {','.join(set(map(str, zones_numbers))) or '[]'}.")
for zone_ndx in range(1,infiltration_timings.shape[1]):
a = wolf2d.inf.array == zone_ndx
b = np.sum(a)
if b == 0:
logging.warning(f"The zone {zone_ndx} doesn't appear in the infiltration map.")
# Compute zones sizes (once so we optimize a bit)
zones_sizes = {}
for zone_ndx in range(1, row.shape[0]):
a = wolf2d.inf.array == zone_ndx
b = np.sum(a)
zones_sizes[zone_ndx] = b
# Spread each infiltration's the water over its zone's meshes.
for row in infiltration_timings:
for zone_ndx in range(1, row.shape[0]):
b = zones_sizes[zone_ndx]
original_inf = row[zone_ndx]
if b > 0:
# Beware it's a surface !
row[zone_ndx] /= b * dxfin * dyfin
#logging.info(f"Infiltration t_start={row[0]:.1f} s, zone {zone_ndx} has {b} cells with {original_inf} m²/s => infiltration per cell = {row[zone_ndx]:.4g} m²/s")
logging.debug("Infiltration timings array after water dispatch")
logging.debug(infiltration_timings)
# print(wolf2d.inf.array)
# print(wolf2d.inf.array.dtype)
# -1 because we must be 0 based (and -xxxx means no zone)
# but wolf is 1 based
# FIXME Align the two
infiltration_zones_array=wolf2d.inf.array.transpose() - 1
# import matplotlib.pyplot as plt
# for row in infiltration_zones_array:
# print(row.tolist())
# plt.imshow(infiltration_zones_array)
# plt.show()
else:
infiltration_timings = None
infiltration_zones_array = None
# ........................................................................
# Admin
# ........................................................................
report_freq = get_report_frequency(wolf2d)
if interpret_freq_type_as_wolf and report_freq.type == SimulationDurationType.SECONDS:
simulation_duration=SimulationDuration.from_seconds(npas)
else:
simulation_duration=SimulationDuration.from_steps(npas)
logging.debug(f"Simulation will run for {npas} steps, records every {int(wolf2d.parameters._writing_frequency)} steps")
if noptpas == 1:
logging.debug("Optimized timestep")
tss = TimeStepStrategy.OPTIMIZED_TIME_STEP
if tstep:
logging.warning("An optimized timestep strategy doesn't need a timestep")
tstep = None
else:
logging.debug("Fixed timestep")
tss = TimeStepStrategy.FIXED_TIME_STEP
# water_height[:,:] = 0
# bathymetry[bathymetry < 0.1] = np.max(bathymetry)
# strong_boundary_cond[:,:] = -99999.0
# strong_boundary_cond[900,:] = 5.0
# strong_boundary_cond[0:200,:] = 0.1
# strong_boundary_cond[1535:1735,:] = 0.1
# strong_boundary_cond[nap == 0] = -99999.0
if glsim is not None:
reload_shaders = froude_max != glsim.froude_max
glsim.set_params(
width=nxfin,
height=nyfin,
dx=dxfin,
dy=dyfin,
ponderation_runge_kutta=ponderation,
courant=vncsouhaite,
froude_max=froude_max,
time_step_strategy=tss,
time_step=tstep,
simulation_duration=simulation_duration,
report_frequency=report_freq,
basefilename=Path(wolf2d.filenamegen).name,
tile_packing_mode=tiles_packing_mode,
optimized_indirection=optimize_indirection,
)
if reload_shaders:
# Have to relaod the shaders if some of their *const*
# have changed... (constants are part of the source code,
# they're not uniform parameters).
logging.debug("Reloading shaders")
glsim._shader_log_path = shader_log_path
else:
glsim = GLSimulation(width=nxfin,
height=nyfin,
dx=dxfin,
dy=dyfin,
ponderation_runge_kutta=ponderation,
courant=vncsouhaite,
froude_max=froude_max,
time_step_strategy=tss,
time_step = tstep,
max_steps=simulation_duration,
report_frequency=report_freq,
basefilename=Path(wolf2d.filenamegen).name,
tile_size=tile_size,
shader_log_path=shader_log_path,
tile_packing_mode=tiles_packing_mode,
optimized_indirection=optimize_indirection)
logging.info("Filling the GPU")
glsim.set_tile_size(tile_size)
glsim.set_buffers(
# Three times the same array of shape (h,w,3)
quantity_arrays=[np.dstack([water_height, qxbin, qybin]) for i in range(3)],
bathymetry_array=bathymetry,
strong_bc_array=strong_boundary_cond, strong_bc_values_array=strong_boundary_cond_values,
manning_array=manning,
active_meshes_array=nap != 0, # FIXME not necessary, remove.
weak_bc_array=weak_boundary_cond,
bc_cells_array=bc_cells_array,
infiltration_zones_array=infiltration_zones_array,
infiltration_timings=infiltration_timings
)
glsim.complete_textures()
# FIXME These info should not be here. Maybe return values ?
if init_gl:
glsim._window_width, glsim._window_height = window_width, window_height
logging.info(f"Model load and GPU data upload took {time.time() - start_time:.1f}s")
return glsim