import sys
if not (
(sys.version_info.major, sys.version_info.minor) == (3, 10) and sys.version_info.micro >= 0
):
print(
f"""I run on python 3.10.x only. You have {sys.version}. Sorry :-(
Please note that Python 3.10 is in security maintenance. So the
latest installable version is 3.10.9. Later 3.10 releases require
you to build python from source yourself."""
)
exit()
from enum import Enum
import logging
from pathlib import Path
from time import time, sleep
from datetime import datetime
from textwrap import wrap
import os
import re
import shutil
import glob
import argparse
import importlib.metadata
import tqdm
import numpy as np
import numpy.ma as ma
import OpenGL
from wolfgpu.tile_packer import TilePackingMode
# Path fiddling to be able to load wolfhece's classes
# if __file__ == "D:\\newgit\\HECEPython\\debug\\test_glsim\\wolfgpu.py":
# _root_dir = Path(__file__).parent.parent.parent
# # Yeah, 1, that's because I want my tests directory to be seen before WolfHECE's one.
# sys.path.insert(1, str(_root_dir))
# else:
# _root_dir = Path(__file__).parent
# sys.path.insert(0, str(_root_dir))
# #print(_root_dir)
# OpenGL.ERROR_CHECKING = False
from OpenGL.GL import GL_BLEND, GL_DEPTH_TEST, glDisable
# from importlib.metadata import version
# wolf_hece_version = tuple(map(int,version('wolfhece').split(".")))
# print(wolf_hece_version[0:2] == (1,8) )
from wolfhece.mesh2d.wolf2dprev import prev_sim2D
# from wolfhece.PyGui import Wolf2DModel # Deprecated in favor of prev_sim2D
from wolfhece.wolf_array import (
WolfArray,
WolfArrayMB,
WolfArrayMNAP,
WOLF_ARRAY_FULL_INTEGER8,
getkeyblock,
WOLF_ARRAY_FULL_SINGLE,
)
from wolfhece.PyVertexvectors import vector, wolfvertex, zone
from wolfgpu.utils import init_global_logging, delete_dir_recursion, nice_timestamp, seconds_to_duration_str
from wolfgpu.test_scenarios import (
save_short_format,
create_steady_flow_swimming_pool_bottom_to_top, scenario1, scenario4,
scenario_cube_drop, scenario_drying, scenario_sine_ground_sine_water,
scenario_small_movement, scenario_still_water, scenario_still_water_rocky_bed
)
from wolfgpu.SimulationRunner import SimulationRunner, MOSTLY_DRY_MESHES_THRESHOLD
from wolfgpu.loaders.wolf_loader import load_sim_to_gpu
from wolfgpu.loaders.simple_sim_loader import load_simple_sim_to_gpu
from wolfgpu.simple_simulation import InfiltrationInterpolation, ReportFrequencyType, SimulationDuration, SimulationDurationType, TimeStepStrategy, SimpleSimulation
from wolfgpu.glsimulation import (
GLSimulation,
plot_evolution_of_simulation,
TileOptimisation,
TIME_SAMPLING_RESOLUTION
)
# from wolfgpu.wolf_utils import force_load # see force_load in 'prev_sim2D' class
# Done here to avoid the message from pygame
import pygame
from wolfgpu.gl_utils import init_pygame, init_glfw
[docs]
VERSION = importlib.metadata.version('wolfgpu')
# resolve(): Make the path absolute, resolving any symlinks.
[docs]
DATA_DIR = Path(__file__).parent.resolve() / "data"
[docs]
WOLF_CLI_EXE = DATA_DIR / "wolfcli.exe"
[docs]
BENCHMARK_AGGREGATED_RESULT_FILE = Path(f"wolfgpu_benchmark_results_{VERSION.replace('.','_')}_{datetime.strftime(datetime.now(),'%Y-%m-%d_%H_%M') }.csv")
# Test case: Crues/2021-07 Vesdre/CSC - Convention - ARNE/Data/Results/Verviers_repository/Simulations/scen_ref_corrected/Q25[50]_not_corctd
# BIG_TEST=Path(__file__).parent / "data" / "Q50_not_corctd"
# Crues\2021-07 Vesdre\CSC - Convention - ARNE\Data\Results\Q25\X1a - Hoegne - Tr 3 - Confluence Wayai - Forges Thiry
[docs]
BIG_TEST = (
DATA_DIR / "X1a - Hoegne - Tr 3 - Confluence Wayai - Forges Thiry"
)
from wolfgpu.wolf_utils import (
write_simulation,
set_report_frequency
)
[docs]
POWER_MODE = (
os.environ["USERNAME"] == "StephaneC" and os.environ["COMPUTERNAME"] == "FEDER1"
)
#POWER_MODE = True
if POWER_MODE:
import cProfile
[docs]
LOGOS = [
r""" __ __ _ __ ___ ___
/ / /\ \ \___ | |/ _| / _ \ / _ \/\ /\
\ \/ \/ / _ \| | |_ / /_\// /_)/ / \ \
\ /\ / (_) | | _/ /_\\/ ___/\ \_/ /
\/ \/ \___/|_|_| \____/\/ \___/""",
r""" ________ __ ___ _______ ______ _______
| | | |.-----.| |.' _| __| __ \ | |
| | | || _ || || _| | | __/ | |
|________||_____||__||__| |_______|___| |_______|""",
r""" __ __ __ _____ __________________ ____ ___
/ \ / \____ | |_/ ____\/ _____/\______ \ | \
\ \/\/ / _ \| |\ __\/ \ ___ | ___/ | /
\ ( <_> ) |_| | \ \_\ \| | | | /
\__/\ / \____/|____/__| \______ /|____| |______/
\/ \/ """,
r""" __ __ _ __ ___ ___ _ _
\ \ / /__| |/ _|/ __| _ \ | | |
\ \/\/ / _ \ | _| (_ | _/ |_| |
\_/\_/\___/_|_| \___|_| \___/"""
]
# def create_model(basefile: str) -> prev_sim2D:
# p = prev_parameters_simul()
# p.nbx = 24
# p.nby = 24
# p.dx = 1.0
# p.dy = 1.0
# model = prev_sim2D(basefile, from_params=p)
# return model
# def show_model_zones(model: prev_sim2D):
# for name, zones in {"model.mymnap.contour": model.mymnap.contour,
# "model.myblocfile.my_vec_blocks": model.bloc_description.my_vec_blocks,
# "model.xyfile.myzones": model.xyfile.myzones,
# "model.mysuxy.myborders": model.sux_suy.myborders,
# "model.myparam.clf.myzones": model.parameters.clf.myzones,
# "model.myparam.clfbx.myzones": model.parameters.clfbx.myzones,
# "model.myparam.clfby.myzones": model.parameters.clfby.myzones,
# }.items():
# print(f"Zones for {name}")
# for zone in zones.myzones:
# print(f" zone: {zone.myname}")
# for polygon in zone.myvectors:
# print(" vector (=polygon)")
# for vertex in polygon.myvertices:
# print(f" x:{vertex.x} y:{vertex.y}")
# def create_sux_suy_for_rectangle(suxy: prev_suxsuy, xmin: int, xmax: int, ymin: int, ymax: int):
# # Coordinates xmin,xmax,ymin,ymax are inclusive
# # SUX/SUY files are used by the VB app
# myparams: prev_parameters_simul
# myparams = suxy.parent.parameters
# dx = myparams.dx
# dy = myparams.dy
# ox = myparams.origx
# oy = myparams.origy
# tx = myparams.translx
# ty = myparams.transly
# suxy.mysux.myvectors.clear()
# suxy.mysuy.myvectors.clear()
# # My hypothesis is that sux borders
# # are set to be on the left border of each
# # cells. So if you want to surround
# # the square:
# # +
# # | y=13
# # +-c-+-d-+
# # f | h y=12
# # +---+---+
# # e | g y=11
# # +-a-+-b-+---+
# # x=1 x=2 x=3
# # You set:
# # SUY = (1,11)a, (2,11)b; (1,13)c, (2,13)d
# # SUX = (1,11)e, (1,12)f; (3,11)g, (3,12)h
# # Horizontal borders (SUY file)
# for y in [ymin, ymax + 1]:
# for k, x in enumerate(range(xmin, xmax + 1)):
# x1 = x * dx + tx + ox
# y1 = y * dy + ty + oy
# vert1 = wolfvertex(x1, y1)
# vert2 = wolfvertex(x1 + dx, y1)
# curborder = vector(name="b" + str(k))
# curborder.add_vertex([vert1, vert2])
# suxy.mysuy.add_vector(curborder)
# # Vertical borders (SUX file)
# for x in [xmin, xmax + 1]:
# for k, y in enumerate(range(ymin, ymax + 1)):
# x1 = x * dx + tx + ox
# y1 = y * dy + ty + oy
# vert1 = wolfvertex(x1, y1)
# vert2 = wolfvertex(x1, y1 + dy)
# curborder = vector(name="b" + str(k))
# curborder.add_vertex([vert1, vert2])
# suxy.mysux.add_vector(curborder)
# def write_sux_suy(suxy: prev_suxsuy, generic_name:Path):
# # --------------------------------------------
# # DEPRECATED use WolfArray.suxsuy_contour
# # --------------------------------------------
# # with open(self.filenamesux, 'r') as f:
# # linesX = f.read().splitlines()
# # with open(self.filenamesuy, 'r') as f:
# # linesY = f.read().splitlines()
# myparams: prev_parameters_simul
# myparams = suxy.parent.myparam
# dx = myparams.dxfin
# dy = myparams.dyfin
# ox = myparams.xminfin
# oy = myparams.yminfin
# tx = myparams.translx
# ty = myparams.transly
# """
# linesX = [np.float64(curline.split(',')) for curline in linesX]
# linesY = [np.float64(curline.split(',')) for curline in linesY]
# k = 1
# for curline in linesX:
# x1 = (curline[0] - 1.) * dx + tx + ox
# y1 = (curline[1] - 1.) * dy + ty + oy
# x2 = x1
# y2 = y1 + dy
# vert1 = wolfvertex(x1, y1)
# vert2 = wolfvertex(x2, y2)
# curborder = vector(name='b' + str(k))
# self.mysux.add_vector(curborder)
# curborder.add_vertex([vert1, vert2])
# k += 1
# """
# with open(generic_name.with_suffix(".sux"), "w") as fout:
# v: vector
# for v in suxy.mysux.myvectors:
# # vert2 doesn't need to be written
# vert1: wolfvertex = v.myvertices[0]
# # vert1 coordinates were computed like this.:
# # x1 = (curline[0] - 1.) * dx + tx + ox
# # y1 = (curline[1] - 1.) * dy + ty + oy
# # We shall invert that:
# # => curline[0] = (x1 -tx -ox)/dx + 1
# x = (vert1.x - tx - ox) / dx + 1
# y = (vert1.y - ty - oy) / dy + 1
# fout.write(f"{x:g},{y:g}\n")
# with open(generic_name.with_suffix(".suy"), "w") as fout:
# for v in suxy.mysuy.myvectors:
# vert1: wolfvertex = v.myvertices[0]
# x = (vert1.x - tx - ox) / dx + 1
# y = (vert1.y - ty - oy) / dy + 1
# fout.write(f"{x:g},{y:g}\n")
# def write_xy(wolf2d: prev_sim2D, generic_path):
# # --------------------------------------------
# # DEPRECATED use WolfArray.suxsuy_contour
# # --------------------------------------------
# from wolfhece.PyVertex import wolfvertex
# t = wolf2d.xyzones.myzones[0].myvectors[0]
# xy = []
# vtx: wolfvertex
# for vtx in t.myvertices:
# xy.append([vtx.x, vtx.y])
# np.savetxt(
# generic_path.with_suffix(".XY"),
# xy,
# header=f"{len(xy):15d}",
# comments="",
# fmt="%12.5f%12.5f",
# )
[docs]
def write_wolf_array(wa, generic_name: Path):
# if's order is important because of inheritance.
if isinstance(wa, WolfArrayMNAP):
# suffix is not in the original filename, setting it ourselves.
wa.filename = str(generic_name.with_suffix(".MNAP"))
wa.write_all()
elif isinstance(wa, WolfArrayMB):
wa: WolfArrayMB
wa.filename = str(generic_name.with_suffix(Path(wa.filename).suffix))
if len(wa.myblocks) >= 1:
wa.write_all()
elif isinstance(wa, WolfArray):
wa.filename = str(generic_name.with_suffix(Path(wa.filename).suffix))
wa.write_all()
# def extend_array(a: WolfArray, x_ext, y_ext):
# # See WolfArray.crop which is the opposite
# assert x_ext >= 0 and y_ext >= 0
# assert isinstance(a, WolfArray)
# assert ma.isMaskedArray(a.array)
# assert a.nbdims == 2, "Only 2D arrays are supported"
# # Remember WolfArrays are masked. Therefore
# # we need to extend mask. In this case, not specifying
# # anything will expand the mask with "dont mask"
# # vaflues.
# # extend vertically
# ex = a.array
# if x_ext > 0:
# # dtype is important: it allows to keep a Fortran friendly
# # type I think.
# ex = ma.append(
# ex,
# np.array([0] * ex.shape[1] * x_ext, dtype=ex.dtype).reshape((x_ext, -1)),
# axis=0,
# )
# a.nbx += x_ext
# # extend horizontally
# if y_ext > 0:
# ex = ma.append(
# ex,
# np.array([0] * ex.shape[0] * y_ext, dtype=ex.dtype).reshape((-1, y_ext)),
# axis=1,
# )
# a.nby += y_ext
# a.array = ex
# a.mask_data(a.nullvalue)
# return a
# def extend_array_mb(a: WolfArrayMB, x_ext, y_ext):
# # See WolfArray.crop which is the opposite
# assert x_ext >= 0 and y_ext >= 0
# assert isinstance(a, WolfArrayMB)
# assert a.nbdims == 2, "Only 2D arrays are supported"
# for block in a.myblocks:
# pass
# a.nbx += x_ext
# a.nby += y_ext
# def make_it_one_block(m: prev_sim2D):
# # ------------------------------------------------------------------------
# # So now we configure the simulation to have (exactly) one block.
# m.init_multibloc()
# # we'll use that block later on to derive all the other arrays
# # The block must be inside the simulation
# # blocki, blockj, block_width, block_height = 8,8,20,10
# # blocki, blockj, block_width, block_height = 1,1,m.myparam.nxfin,m.myparam.nyfin
# blocki, blockj, block_width, block_height = (
# 3,
# 3,
# m.myparam.nxfin - 4,
# m.myparam.nyfin - 4,
# )
# # It seems that block numbering starts at 1
# # (see write_txt_header in WolfArray writer)
# block_number = 1
# block_key = getkeyblock(block_number, addone=False)
# first_block_array = WolfArray(whichtype=WOLF_ARRAY_FULL_INTEGER8)
# # FIXME idx should be set on the cosntructor and should
# # be immutable ('cos it is there to keep structure)
# first_block_array.idx = block_key
# first_block_array.nbx = block_width
# first_block_array.nby = block_height
# first_block_array.dx = m.parameters.dx
# first_block_array.dy = m.parameters.dy
# first_block_array.origx = blocki # Relative to the MNAP origin
# first_block_array.origy = blockj # Relative to the MNAP origin
# first_block_array.translx = m.myparam.translx
# first_block_array.transly = m.myparam.transly
# first_block_array.isblock = True
# first_block_array.blockindex = 0
# ones = np.ones((block_width, block_height), order="F", dtype=np.float32)
# first_block_array.array = ma.MaskedArray(ones, mask=np.zeros_like(ones))
# # MNAP array has two sides:
# # - it has an array the size of the whole computation domain
# # which holds nieghbours references (???) of blocks (self.array)
# # - it has an array that represents the NAP of each block (self.myblocks)
# # FIXME Is that even correct ?
# zeros = np.zeros((m.myparam.nxfin, m.myparam.nyfin), order="F", dtype=np.float32)
# m.mymnap.array = ma.array(zeros, mask=np.zeros_like(zeros))
# m.mymnap.array[0:5, :] = 1
# m.mymnap.origx = 0.0
# m.mymnap.origy = 0.0
# m.mymnap.add_block(first_block_array)
# # For MNAP array there's an extra step for catching the contour
# # of the block
# block_bounding_rectangle = vector(name=f"{block_key}-BoundingBox")
# block_bounding_rectangle.add_vertex(
# wolfvertex(*m.mymnap.get_xy_from_ij(i=0, j=0, which_block=block_number))
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(i=m.myparam.nxfin, j=0, which_block=block_number)
# )
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(
# i=m.myparam.nxfin, j=m.myparam.nyfin, which_block=block_number
# )
# )
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(i=0, j=m.myparam.nyfin, which_block=block_number)
# )
# )
# # FIXME Don't use the block key as a name as sooner or later
# # someone will use that name to index dictionaries...
# bloc_zone = zone(name=block_key)
# bloc_zone.add_vector(block_bounding_rectangle)
# block_bounding_rectangle.parentzone = bloc_zone
# m.mymnap.contour.add_zone(bloc_zone)
# # ------------------------------------------------------------------------
# # Create the various block's arrays
# def make_simple_block_like(mold: WolfArray, key):
# # with mold, we copy dimensions and contents.
# bloc = WolfArray(whichtype=WOLF_ARRAY_FULL_SINGLE, mold=mold, preload=False)
# # FIXME mold copies datra but since our mold has no
# # array, it fails and that's what I fix here:
# bloc.array = ma.MaskedArray(
# # np.zeros((bloc.nbx, bloc.nby), order='F', dtype=np.float32))
# np.zeros(
# (m.myparam.nxfin + 4, m.myparam.nxfin + 4), order="F", dtype=np.float32
# )
# )
# # FIXME idx should be renamed "key" since it is built
# # by the function `getkeyblock`
# bloc.idx = key
# bloc.mask_reset() # Unmask everything
# return bloc
# m.topini.add_block(make_simple_block_like(first_block_array, block_key))
# m.topini.head_blocks[block_key] = first_block_array.get_header()
# m.hbinb.add_block(make_simple_block_like(first_block_array, block_key))
# m.hbinb.head_blocks[block_key] = first_block_array.get_header()
# m.qxbinb.add_block(make_simple_block_like(first_block_array, block_key))
# m.qxbinb.head_blocks[block_key] = first_block_array.get_header()
# m.qybinb.add_block(make_simple_block_like(first_block_array, block_key))
# m.qybinb.head_blocks[block_key] = first_block_array.get_header()
# # ------------------------------------------------------------------------
# general: zone = m.myblocfile.my_vec_blocks[0]
# # Define the polygon bounding the simulation. The polygon has
# # whaterver shape, but we choose a rectangle because it's easier.
# # Be aware of the fact that coordinates conversion sets the border in the
# # middle of the cells.
# # FIXME Typechecker fails: class bloc_external_border(vector):
# bounding_rectangle: block_contour = general[0]
# if bounding_rectangle.nbvertices == 0:
# bounding_rectangle.add_vertex(
# wolfvertex(*m.mymnap.get_xy_from_ij(i=0, j=0, which_block=block_number))
# )
# bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(
# i=m.myparam.nxfin, j=0, which_block=block_number
# )
# )
# )
# bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(
# i=m.myparam.nxfin, j=m.myparam.nyfin, which_block=block_number
# )
# )
# )
# bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(
# i=0, j=m.myparam.nyfin, which_block=block_number
# )
# )
# )
# bounding_rectangle.close_force()
# bounding_rectangle.set_limits()
# else:
# # Already initialized somewhere else
# pass
# # Now the blocks, one by one
# # The block bounding rectangle will be shared.
# block_bounding_rectangle = block_contour(
# name=f"{block_key}-ExternalBorder", lines=None
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(*m.mymnap.get_xy_from_ij(i=0, j=0, which_block=block_number))
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(i=m.myparam.nxfin, j=0, which_block=block_number)
# )
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(
# i=m.myparam.nxfin, j=m.myparam.nyfin, which_block=block_number
# )
# )
# )
# block_bounding_rectangle.add_vertex(
# wolfvertex(
# *m.mymnap.get_xy_from_ij(i=0, j=m.myparam.nyfin, which_block=block_number)
# )
# )
# block_bounding_rectangle.close_force()
# block_bounding_rectangle.set_limits()
# bfile = blocks_file(parent=m)
# bextent = block_description(parent=bfile, idx=0, lines=None)
# # Set the magnetic grid min/max
# # FIXME Is this correct ?
# bextent.nb = block_bounding_rectangle.nbvertices
# bextent.xmin = block_bounding_rectangle.minx
# bextent.xmax = block_bounding_rectangle.maxx
# bextent.ymin = block_bounding_rectangle.miny
# bextent.ymax = block_bounding_rectangle.maxy
# bfile.my_blocks.append(bextent)
# all_block_extents: zone = m.myblocfile.my_vec_blocks[1]
# all_block_extents.add_vector(block_bounding_rectangle)
# # ------------------------------------------------------------------------
# for bext in m.myblocfile.my_blocks:
# bext: block_description
# b = bext.external_border
# # mini, minj = m.napbin.get_ij_from_xy(b.minx, b.miny)
# # maxi, maxj = m.napbin.get_ij_from_xy(b.maxx, b.maxy)
# for i in range(block_width):
# for j in range(block_height):
# m.mymnap.array[i, j] = 1
# write_simulation(m)
# def set_block_extent(wolf2d: prev_sim2D, polygon: list):
# # Polygon is given in matrix coordinates, with (1,1) origin.
# # (so it's not real coordiantes, it's no float and it doesn't
# # start at (0,0))
# assert len(polygon) >= 3, "Polygon must have at least 3 vertices"
# # Close the loop
# polygon = [x for x in polygon]
# polygon.append(polygon[0])
# """
# The .bloc file is made like this:
# See class blocks_file method write_file.
# """
# # Block extents come in several places:
# # See class bloc_external_border and how it is used by class blocks_file.
# # - Zone 0 is the general zone, with one arbitrary enclosing polygon
# # - Zone 1 is the list of blocks extent (rectangular polygon), this
# # one is accompanied by the blocks_file.my_blocks list
# # model.my_blocks
# # [i] + block_extent() objects (1 per block)
# # model.blocks_file.my_vec_blocks.myzones (class: zone)
# # [0] + General Zone (here called: general_zone)
# # + external border (.myvectors[0])
# # [1] + Block extents (here called: blocks_extent_zone)
# # +-> myvectors.append(new_bloc_extent.external_border)
# # +-> myvectors.append(block_ext_border)
# general_zone: zone = wolf2d.myblocfile.my_vec_blocks.myzones[0]
# blocks_extent_zone: zone = wolf2d.myblocfile.my_vec_blocks.myzones[1]
# assert (
# len(general_zone.myvectors) == 1
# ), "Hypothesis: There's only one contour in the general zone."
# assert isinstance(general_zone.myvectors[0], block_contour)
# # .......................................................
# general_zone.myvectors.clear()
# # Although bloc_external_border() can be built with a list of text lines,
# # it's not done like that. It is done in blocks_file class because the
# # syntax of the lines is different (negative vertex count; 3 columns, etc.)
# assert wolf2d.myblocfile.interior == True
# # FIXME vertices (parent of bloc_external_border) needs
# # a name to be able to receive new vertex. If no name is
# # provided, then vector.vertices[] array is not init'd
# # and one cannot append anything to it.
# block_ext_border = block_contour(is2D=True, name="external border")
# trlx = wolf2d.myparam.translx
# trly = wolf2d.myparam.transly
# dx = wolf2d.parameters.dx
# dy = wolf2d.parameters.dy
# # When read in blocks_file.read_file(), the coordinates are
# # added with trlx/trly BUT not multiplied by dx/dy
# # cx, cy are given in cell numbers : [1,...]
# # we translate to real coordinates
# polygon_crd = [
# (float(cx - 1) * dx + trlx, float(cy - 1) * dy + trly) for cx, cy in polygon
# ]
# # This is the black one
# for cx, cy in polygon_crd:
# # x, y = wolf2d.napbin.get_xy_from_ij(cx, cy, aswolf=True)
# x, y = cx, cy
# curvert = wolfvertex(x, y, 0.0) # Doesn't xform anything
# block_ext_border.add_vertex(curvert)
# general_zone.add_vector(block_ext_border)
# general_zone.find_minmax(update=True)
# # Now the block extents ............................................
# assert (
# len(blocks_extent_zone.myvectors) == wolf2d.myparam.nblocs
# ), "an hypothesis I keep in mind, still not sure it's correct though"
# assert len(wolf2d.myblocfile.my_blocks) == wolf2d.myparam.nblocs
# # assert isinstance(blocks_extent_zone.myvectors[0], bloc_external_border)
# # assert isinstance(wolf2d.myblocfile.my_blocks[0], block_description)
# # This is only one block
# # FIXME ? bloc_exten() reads: x + tx (== parent.parent.myparam.translx)
# # This is the red one
# t = [str(len(polygon))] + [f"{x-1},{y-1}" for x, y in polygon]
# ##emprise du bloc en accord avec le grid magnétique
# a = np.array(polygon) - 1
# t.append(f"{np.min(a[:,0])},{np.max(a[:,0])}") # xmin, xmax
# t.append(f"{np.min(a[:,1])},{np.max(a[:,1])}") # ymin, ymax
# # This call will set up the .external_border member according
# # to the lines parameter
# new_bloc_extent = block_description(parent=wolf2d.myblocfile, lines=t)
# # FIXME Fishy! parentzone is declares as "Zones" but it is
# # sometimes assigned with "Zone" (singular)...
# new_bloc_extent.external_border.parentzone = blocks_extent_zone
# # FIXME This is rather strange, block_description discards param's dx,dy !!!
# # Maybe it's because it counts in grid units ?
# new_bloc_extent.dx = wolf2d.parameters.dx
# new_bloc_extent.dy = wolf2d.parameters.dy
# # ......................................................
# # Updating my_blocks with new_bloc_extent
# wolf2d.myblocfile.my_blocks.clear()
# wolf2d.myblocfile.my_blocks.append(new_bloc_extent)
# # ......................................................
# # Updating blocks_extent_zone.myvectors with new_bloc_extent
# blocks_extent_zone.myvectors.clear()
# blocks_extent_zone.myvectors.append(new_bloc_extent.external_border)
# # ......................................................
# # Updating blocks_extent_zone with bloc_external_border
# # Although bloc_external_border() can be built with a list of text lines,
# # it's not done like that. It is done in blocks_file class because the
# # syntax of the text lines is different (negative vertex count; 3 columns, etc.)
# # FIXME vertices (parent of bloc_exteranl_border) needs
# # a name to be able to receive new vertex. If no name is
# # provided, then vector.vertices[] array is not init'd
# # and one cannot append anything to it.
# # FIXME ??? When this one reads its input (if provided),
# # it makes no conversion on vertices coordinates.
# block_ext_border = block_contour(is2D=True, name="external border")
# for x, y in polygon:
# curvert = wolfvertex(x, y, 0.0)
# block_ext_border.add_vertex(curvert)
# blocks_extent_zone.add_vector(block_ext_border)
# blocks_extent_zone.find_minmax(update=True)
# def extend_mesh(wolf2d: prev_sim2D, x_ext: float, y_ext: float):
# # x_ext, y_ext: number (integer) of meshes to add
# if x_ext == 0 and y_ext == 0:
# return
# print(
# f"{wolf2d.napbin.nbx}x{wolf2d.napbin.nby} numpy:{wolf2d.napbin.array.shape} !null:{wolf2d.napbin.nbnotnull}"
# )
# extend_array(wolf2d.napbin, x_ext, y_ext)
# extend_array(wolf2d.top, x_ext, y_ext)
# extend_array(wolf2d.hbin, x_ext, y_ext)
# extend_array(wolf2d.frot, x_ext, y_ext)
# extend_array(wolf2d.qxbin, x_ext, y_ext)
# extend_array(wolf2d.qybin, x_ext, y_ext)
# extend_array_mb(wolf2d.hbinb, x_ext, y_ext)
# extend_array_mb(wolf2d.qxbinb, x_ext, y_ext)
# extend_array_mb(wolf2d.qybinb, x_ext, y_ext)
# extend_array_mb(wolf2d.topini, x_ext, y_ext)
# wolf2d.myparam.nxfin += x_ext
# wolf2d.myparam.nyfin += y_ext
# print(
# f"{wolf2d.napbin.nbx}x{wolf2d.napbin.nby} numpy:{wolf2d.napbin.array.shape} !null:{wolf2d.napbin.nbnotnull}"
# )
[docs]
class Scenario(Enum):
[docs]
SWIMMING_POOL_WITH_HOLE = 2
[docs]
def bbox2(img):
# From https://stackoverflow.com/questions/31400769/bounding-box-of-numpy-array
rows = np.any(img, axis=1)
cols = np.any(img, axis=0)
rmin, rmax = np.where(rows)[0][[0, -1]]
cmin, cmax = np.where(cols)[0][[0, -1]]
return rmin, rmax, cmin, cmax
[docs]
def bounding_polygon_xy(wa: WolfArray):
# Bounding polygon is for *borders* which are
# located on the left or bottom of each cell.
rmin, rmax, cmin, cmax = bbox2(wa.array)
# The coord conversion adds a 0.5 when covnerting. So
# we adapt our coordinates so we have a bouding polygon
# which is outside the data it surrounds.
polygon = wa.get_xy_from_ij_array(
np.array([(cmin - 1, rmin - 1), (cmin, rmax), (cmax, rmax), (cmax, rmin - 1)])
)
return polygon
[docs]
def mask_border(wa: WolfArray):
""" Set a new mask applying only to a one cell around
the matrix. It is usefule because SUXY is not smart
enough to handle borders which are on the very border
of matrices."""
wa.mask_reset()
m = np.zeros_like(wa.array, dtype=bool)
m[0, :] = True
m[:, 0] = True
m[m.shape[0] - 1, :] = True
m[:, m.shape[1] - 1] = True
wa.array.mask = m
wa.nbnotnull = wa.array.count()
# def copy_simulation(
# source_path,
# destination_path,
# generic_name,
# scenario=Scenario.TWO_ROOMS,
# scenario_param=None,
# niter=1000,
# ) -> prev_sim2D:
# start_time = time()
# logging.info(f"Copying model data form {source_path}")
# generic_path = destination_path / generic_name
# wolf2d: prev_sim2D = prev_sim2D(source_path, splash=False, in_gui=False)
# force_load(wolf2d)
# # wolf2d = create_model(str(generic_path))
# show_model_zones(wolf2d)
# wolf2d.myparam.npas = niter
# if scenario == Scenario.LINEAR:
# wolf2d.myparam.ntypefreq = 0 # en pas ?
# wolf2d.myparam.noptpas = 0 # =1 si optimisation du pas de temps
# wolf2d.myparam.dur = 0.01 # durée pas de temps
# wolf2d.myparam.ntyplimit = (
# 0 # 0 si pas de limiteur, 1 si barth jesperson, 2 si venkatakrishnan
# )
# wolf2d.myparam.my_param_blocks[
# 0
# ].nconflit = 1 ##< gestion des conflits (0), ou juste en centré (1) ou pas (2)
# wolf2d.myparam.translx = 0.0
# wolf2d.myparam.transly = 0.0
# wolf2d.myparam.xminfin = 0.0
# wolf2d.myparam.yminfin = 0.0
# wolf2d.napbin.origx = wolf2d.myparam.xminfin
# wolf2d.napbin.origy = wolf2d.myparam.yminfin
# wolf2d.napbin.translx = wolf2d.myparam.translx
# wolf2d.napbin.transly = wolf2d.myparam.transly
# # Pattern used to align WolfCLI and GPUSim
# pattern = [
# 0.5,
# 1,
# 0.5,
# 1,
# 1,
# 1,
# 0,
# 0,
# 1,
# 1,
# 0,
# 0,
# 1,
# 1,
# 0,
# 1,
# 0,
# 1,
# ] # ,1,0.1,1]
# pattern = [10] * len(pattern)
# # Putting 0 height blows some masks in the data model.
# pattern = [max(0.2, x) for x in pattern]
# # xstart, xend coordinates are in grid coord (one based)
# xstart = 4 # 1-based
# xend = xstart + len(pattern) - 1 # xend is inclusive.
# # y's covered by the swimming pool, in mesh coordinates (1-based)
# ys = [(wolf2d.myparam.nyfin // 2) + i for i in range(-5, +5 + 1)]
# ystart, yend = ys[0], ys[-1]
# polygon = [
# (xstart - 3, ystart - 2),
# (xstart - 3, yend + 2),
# (xend + 3, yend + 2),
# (xend + 3, ystart - 2),
# ]
# set_block_extent(wolf2d, polygon)
# general_zone: zone = wolf2d.myblocfile.my_vec_blocks.myzones[0]
# blocks_extent_zone: zone = wolf2d.myblocfile.my_vec_blocks.myzones[1]
# vec_poly = vector(name="clip")
# for cx, cy in polygon:
# x, y = wolf2d.napbin.get_xy_from_ij(cx, cy, aswolf=True)
# vec_poly.add_vertex(wolfvertex(x, y, 0.0))
# wolf2d.napbin.mask_reset()
# wolf2d.napbin.array[:, :] = -1
# wolf2d.napbin.mask_outsidepoly(vec_poly)
# # wolf2d.top.mask_outsidepoly(general_zone.myvectors[0])
# # wolf2d.hbin.mask_outsidepoly(general_zone.myvectors[0])
# # wolf2d.qxbin.mask_outsidepoly(general_zone.myvectors[0])
# # wolf2d.qybin.mask_outsidepoly(general_zone.myvectors[0])
# # wolf2d.frot.mask_outsidepoly(general_zone.myvectors[0])
# # ....................................................................
# # Lets fill the topography with some concrete, 10m high.
# wolf2d.top.mask_data(-1)
# wolf2d.top.array[:, :] = 10
# # Now we dig a swimming pool along the x axis.
# for y in ys:
# # Here we access an array, so x,y are zero based (intead of one based)
# # xend gets a +1 to be inclusive
# wolf2d.top.array[xstart - 1 : xend - 1 + 1, y - 1] = 0
# # ....................................................................
# # Set initial water height (.hbin file)
# # wolf2d.hbin.mask_data(-1) # mask nothing
# wolf2d.hbin.reset()
# # Put water in the swimming pool
# for y in ys:
# wolf2d.hbin.array[xstart - 1 : xend - 1 + 1, y - 1] = pattern
# # ....................................................................
# # Generate weak boundary conditions
# # Those will appear in the .PAR file
# xbc = xstart
# if True:
# ys_of_bc = [ys[len(ys) // 2]]
# wolf2d.myparam.impfbxgen = len(ys_of_bc) # Seems not used in PAR file ???
# for y in ys_of_bc:
# # Fixing coordinates so that everything
# # is aligned and VB is happy...
# # Fortran fix adds +2,+1, here we undo this.
# i, j = xbc - 2, y - 1
# wolf2d.myparam.clfbx.mybc.append(
# prev_boundary_condition(i, j, BoundaryConditionsTypes.QX.value, 1.0)
# )
# # ....................................................................
# # Generate strong boundary conditions
# # FIXME ??? These seem unsupported by the python data model right now
# # FIXME Moreover, this prevents WolfCLI from running.
# if False:
# wolf2d.myparam.impfgen = len(ys)
# for y in ys:
# wolf2d.myparam.clf.mybc.append(
# prev_boundary_condition(
# xbc, y, BoundaryConditionsTypes.QX.value, 1.0
# )
# ) # x,y mesh; indicator (2 == qx), value
# # ....................................................................
# # hbinb seems empty
# # (FIXME because there's only one block ? and binb means bin multiblocks)
# assert wolf2d.hbinb.array == None
# assert len(wolf2d.hbinb.myblocks) == 0
# # ....................................................................
# assert wolf2d.topini.array == None
# assert len(wolf2d.topini.myblocks) == 0
# # assert wolf2d.topinifine.array == None
# # Polygon structure is:
# # 2---3
# # | |
# # 1---4
# xstart, xend = xstart, xend
# ystart, yend = ys[0], ys[-1]
# # create_sux_suy_for_rectangle(wolf2d.mysuxy, xstart, xend, ystart, yend)
# # FIXME These are mysterious corrections
# # Without these, WolfCli doesn't align its grids
# # onto ours.
# # Keep in mind these may depends on translations
# # If one opens a sim in the VB app, one cans see
# # that these translate to a 2-cells border around the
# # matrix area of interest.
# # I figured out all thes values while setting
# # all (to the best of my knowledge) translations to zero.
# xstart -= 2
# xend += 2
# ystart -= 2
# yend += 2
# elif scenario == Scenario.TWO_ROOMS:
# x_ext, y_ext = scenario_param or 100, scenario_param or 100
# new_x, new_y = wolf2d.myparam.nxfin + x_ext, wolf2d.myparam.nyfin + y_ext
# corr = 6
# polygon = [
# (1, 1),
# (1, wolf2d.myparam.nyfin + y_ext - corr),
# (wolf2d.myparam.nxfin + x_ext - corr, wolf2d.myparam.nyfin + y_ext - corr),
# (wolf2d.myparam.nxfin + x_ext - corr, 1),
# ]
# set_block_extent(wolf2d, polygon)
# extend_mesh(wolf2d, x_ext, y_ext)
# xmid = new_x // 2
# wolf2d.top.mask_data(-1)
# wolf2d.top.array[:, :] = 10
# wolf2d.top.array[4 : (new_x - 4), 4 : (new_y - 4)] = 0.1
# # Wall cuttin the pool in two
# wolf2d.top.array[xmid, :] = 10
# # Little hole
# wolf2d.top.array[xmid, new_y // 2 : new_y // 2 + 2] = 0.1
# wolf2d.hbin.mask_data(-1)
# wolf2d.hbin.array[:, :] = 0.1
# wolf2d.hbin.array[0:xmid, :] = 1
# wolf2d.hbin.array[:, new_y - 10] = 1
# wolf2d.hbin.array[new_x - 10, :] = 1
# wolf2d.hbin.array[wolf2d.top.array > 5] = 0
# wolf2d.myparam.npas = 200
# wolf2d.myparam.ntypefreq = 0 # en pas ?
# wolf2d.myparam.noptpas = 0 # =1 si optimisation du pas de temps
# wolf2d.myparam.dur = 0.01 # durée pas de temps
# wolf2d.myparam.ntyplimit = (
# 0 # 0 si pas de limiteur, 1 si barth jesperson, 2 si venkatakrishnan
# )
# wolf2d.myparam.my_param_blocks[
# 0
# ].nconflit = 1 ##< gestion des conflits (0), ou juste en centré (1) ou pas (2)
# elif scenario == Scenario.TWO_SMALL_ROOMS:
# x_ext, y_ext = 1, 1
# new_x, new_y = wolf2d.myparam.nxfin + x_ext, wolf2d.myparam.nyfin + y_ext
# corr = 2 * 3 - 1
# polygon = [
# (1, 1),
# (1, wolf2d.myparam.nyfin + y_ext - corr),
# (wolf2d.myparam.nxfin + x_ext - corr, wolf2d.myparam.nyfin + y_ext - corr),
# (wolf2d.myparam.nxfin + x_ext - corr, 1),
# ]
# set_block_extent(wolf2d, polygon)
# extend_mesh(wolf2d, x_ext, y_ext)
# width, height = wolf2d.myparam.nxfin, wolf2d.myparam.nyfin
# xmid = width // 2
# wolf2d.top.mask_data(-1)
# wolf2d.top.array[:, :] = 10
# wolf2d.top.array[4 : (width - 4), 4 : (height - 4)] = 0.1
# # Wall cuttin the pool in two
# wolf2d.top.array[xmid, :] = 10
# # Little hole
# wolf2d.top.array[xmid, height // 2 : height // 2 + 2] = 0.1
# wolf2d.hbin.mask_data(-1)
# wolf2d.hbin.array[:, :] = 0.5
# wolf2d.hbin.array[0:xmid, :] = 1
# wolf2d.hbin.array[wolf2d.top.array > 5] = 0
# wolf2d.myparam.ntypefreq = 0 # en pas ?
# wolf2d.myparam.noptpas = 0 # =1 si optimisation du pas de temps
# wolf2d.myparam.dur = 0.005 # durée pas de temps
# wolf2d.myparam.ntyplimit = (
# 0 # 0 si pas de limiteur, 1 si barth jesperson, 2 si venkatakrishnan
# )
# wolf2d.myparam.my_param_blocks[
# 0
# ].nconflit = 1 ##< gestion des conflits (0), ou juste en centré (1) ou pas (2)
# elif scenario == Scenario.SWIMMING_POOL_WITH_HOLE:
# wolf2d.top.mask_data(-1)
# wolf2d.top.array[:, :] = 10
# wolf2d.top.array[4 : (new_x - 4), 4 : (new_y - 4)] = 0.1
# wolf2d.hbin.mask_data(-1)
# wolf2d.hbin.array[:, :] = 0.1
# cx = new_x // 2
# cy = new_y // 2
# dx = new_x // 5
# dy = new_y // 5
# wolf2d.hbin.array[cx - dx : cx + dx, cy - dy : cy + dy] = 1
# wolf2d.hbin.array[wolf2d.top.array > 5] = 0
# wolf2d.myparam.npas = 10000
# wolf2d.myparam.ntypefreq = 0 # en pas ?
# wolf2d.myparam.noptpas = 0 # =1 si optimisation du pas de temps
# wolf2d.myparam.dur = 0.01 # durée pas de temps
# wolf2d.myparam.ntyplimit = (
# 0 # 0 si pas de limiteur, 1 si barth jesperson, 2 si venkatakrishnan
# )
# wolf2d.myparam.my_param_blocks[
# 0
# ].nconflit = 1 ##< gestion des conflits (0), ou juste en centré (1) ou pas (2)
# wolf2d.frot.array[:, :] = 0.2
# write_simulation(wolf2d, source_path, generic_path, destination_path, generic_name)
# # print("-"*80)
# # print(wolf2d.top.array)
# # print("-"*80)
# # print(wolf2d.hbin.array)
# # print("-"*80)
# # print(wolf2d.hbin.array + wolf2d.top.array)
# # print("-"*80 + "Frottement")
# # print(wolf2d.frot.array)
# # print(wolf2d.myparam.to_yaml())
# return wolf2d
[docs]
def run_wolf_cli(destination_dir: Path):
# destination_dir: wher WolfCLI is located
# Return a handler on which one can .wait()
import subprocess
import shlex
wcp = str(destination_dir / "wolfcli.exe")
if os.path.exists(wcp):
cmd = wcp + " run_wolf2d_prev genfile=simul"
logging.info(
f"Running wolf cli in {destination_dir}, as {' '.join(shlex.split(cmd, posix=False))}"
)
return subprocess.Popen(
shlex.split(cmd, posix=False),
cwd=destination_dir,
stdout=subprocess.DEVNULL,
)
else:
logging.error(
f"Can't run wolfcli.exe because I can't find it at {destination_dir}"
)
return None
[docs]
def force_cli_params_in_wolf_model(cli_args: argparse.Namespace, m: prev_sim2D):
if cli_args.freq is not None:
try:
set_report_frequency(m, cli_args.freq)
except ValueError as ex:
logging.warning(f"I can't set the reporting frequency {cli_args.freq} into a Wolf model. This will be handled only on the GPU side.")
if cli_args.npas is not None:
m.parameters._nb_timesteps = int(cli_args.npas)
if m.parameters._writing_mode == 1:
logging.warning("You gave -npas but the work Wolf scenario has type frequency report = seconds. -npas will be interpreted as seconds by WolfCLI.")
elif cli_args.sim_duration is not None:
m.parameters._nb_timesteps = SimulationDuration.from_seconds(cli_args.sim_duration).duration
if m.parameters._writing_mode == 0:
logging.warning("You gave -sim-duration but the Wolf scenario has type frequency report = iterations. -sim-duration will interpreted as a number of iterations by WolfCLI.")
if cli_args.optim_pas is not None:
m.parameters._scheme_optimize_timestep = int(cli_args.optim_pas)
if cli_args.dur is not None:
m.parameters._timestep_duration = float(cli_args.dur)
p = cli_args.ponderation
if p is not None:
if p == "euler":
m.parameters._scheme_rk = 1.
# logging.warning(
# "I don't know how to set Euler-only on a Wolf simulation. But I do know on a GPU simulation, so I'll do that"
# )
else:
m.parameters._scheme_rk = p
if cli_args.courant is not None:
m.parameters._scheme_cfl = cli_args.courant
if cli_args.manning is not None:
m.frot.array[:, :] = cli_args.manning
if cli_args.froude_max:
assert cli_args.froude_max > 0.0, "Froude limit should be > 0"
m.parameters.blocks[0]._froude_max = cli_args.froude_max
[docs]
def force_params_in_gpu(cli_args: argparse.Namespace, gpu: GLSimulation):
# Change parameters which can only be set in the GPU model (not in the Wolf model)
# Well, we start with an exception: I set wolf parameters here. Why ?
# Because those parameters are more flexible in the GPU than in Wolf.
if cli_args.dur:
gpu.set_fixed_time_step( cli_args.dur)
if cli_args.froude_max:
gpu.set_froude_limit(cli_args.froude_max)
if cli_args.courant:
gpu.set_courant( cli_args.courant)
if cli_args.freq is not None:
gpu.set_report_frequency( SimulationDuration.from_str(cli_args.freq))
if cli_args.npas is not None:
gpu.set_simulation_duration( SimulationDuration.from_steps(cli_args.npas))
elif cli_args.sim_duration is not None:
gpu.set_simulation_duration( SimulationDuration.from_seconds(cli_args.sim_duration))
if cli_args.optim_pas is not None:
if int(cli_args.optim_pas) == 0:
gpu.time_step_strategy = TimeStepStrategy.FIXED_TIME_STEP
else:
gpu.time_step_strategy = TimeStepStrategy.OPTIMIZED_TIME_STEP
# if cli_args.freq is not None:
# gpu.report_frequency = cli_args.freq
p = cli_args.ponderation
if p is not None:
if p == "euler":
gpu.set_euler_ponderation()
else:
gpu.set_rk_ponderation(p)
if cli_args.infilerp is None:
gpu.set_infiltration_interpolation(InfiltrationInterpolation.NONE)
elif cli_args.infilerp == "linear":
gpu.set_infiltration_interpolation(InfiltrationInterpolation.LINEAR)
else:
raise Exception(f"Unrecognized infiltration interpolation strategy: '{cli_args.infilerp}'")
if cli_args.dry_up_loops is not None:
assert cli_args.dry_up_loops >= 0, "If you specify the number of dry up loops then it must be >= 0"
gpu.set_total_number_of_dry_up_outer_loops(cli_args.dry_up_loops)
if cli_args.sim_duration:
gpu.set_simulation_duration(SimulationDuration.from_seconds(cli_args.sim_duration))
[docs]
def print_sim_params(sim: GLSimulation, runner: SimulationRunner,
manning: np.array, active: np.array, h: np.array):
mask = active != 0
total_wet = np.count_nonzero(h[mask] > 0)
m = manning[mask]
total_active = np.sum(mask)
domain = sim.width*sim.height
print("\nSimulation parameters")
print("---------------------\n")
print(f"Start time : {nice_timestamp()}")
print(f"Domain : {sim.width}x{sim.height} meshes (dx={sim.dx:g} m.,dy={sim.dx:g} m.), {100*(total_active/domain):.1f}% active, {100*(total_wet/domain):.1f}% wet")
if sim.time_step_strategy == TimeStepStrategy.FIXED_TIME_STEP:
print(f"Timestep : Fixed, {sim.time_step:g} s.")
else:
print(f"Timestep : Optimized")
if sim.report_frequency is not None:
rep = f"1 report every {sim.report_frequency}"
if sim.report_frequency.type == SimulationDurationType.SECONDS:
rep += f" (time sampled every {TIME_SAMPLING_RESOLUTION} steps)"
else:
rep = "no reports recorded"
print(f"Duration : {sim.simulation_duration}, {rep}")
if sim.ponderation_runge_kutta:
print(
f"Ponderation : Runge-Kutta {sim.ponderation_runge_kutta:g}/{1-sim.ponderation_runge_kutta:g}"
)
else:
print("Ponderation : Euler")
if sim.infiltration_interpolation == InfiltrationInterpolation.NONE:
interp = "None"
elif sim.infiltration_interpolation == InfiltrationInterpolation.LINEAR:
interp = "Linear"
print(f"Infiltration interpolation : {interp}")
# FIXME sim._tile_size: abstraction leak !
print(f"Tile size : {sim._tile_size}x{sim._tile_size} meshes")
print(f"Courant number : {sim.courant:g}") # FIXME Non dimensional, correct ?
print(f"Froude limit : {sim.froude_max:g}")
print(f"g : {sim.g:g} m/s²")
if sim.total_number_of_dry_up_outer_loops is None:
dry_up = "Solve until done"
elif sim.total_number_of_dry_up_outer_loops == 0:
dry_up = "Disabled"
elif sim.total_number_of_dry_up_outer_loops == 1:
dry_up = "One iteration"
else:
dry_up = f"{sim.total_number_of_dry_up_outer_loops} iterations"
print(f"Dry up handling: {dry_up}")
print(f"Tiles packing : {sim._tile_packing_mode}; {['standard','optimized'][sim.optimize_tile_indirection]} tile indirection")
print(f"Mean of Manning coefficients: {np.mean(m):.3f} [s.m^{{-1/3}}]")
print(f"Early out threshold: {runner.early_out_threshold}")
print(f"Mostly dry mesh threshold: {MOSTLY_DRY_MESHES_THRESHOLD} m")
# #print(active)
# print(np.histogram(m))
# from matplotlib import pyplot as plt
# #plt.hist(m)
# plt.plot( np.sort(m))
# plt.show()
# #print(m)
print()
[docs]
def benchmark(destination_dir, cli_args):
from wolfgpu.SimulationRunner import estimate_best_window_size
if False and POWER_MODE:
SIM_SIZES = [512,1024,2048]
else:
SIM_SIZES = [512,512,512,1024,1024,1024,2048,2048,2048,1024*4,1024*5,1024*6,1024*7,1024*8]
if cli_args.ap is None:
cli_args.ap = 0.8
# Make sure we measure performance of the simulator and not
# of the reporting.
cli_args.freq = "100000s"
from OpenGL.GL import glGetIntegerv, GL_MAJOR_VERSION, GL_MINOR_VERSION, glGetString, GL_VERSION, GL_VENDOR, GL_SHADING_LANGUAGE_VERSION, GL_RENDERER
with open(BENCHMARK_AGGREGATED_RESULT_FILE,"w") as res_out:
pygame.init()
pygame.display.set_mode((100, 100), pygame.OPENGL|pygame.DOUBLEBUF)
glnfo = f"\"OpenGl version: {glGetIntegerv(GL_MAJOR_VERSION)}.{glGetIntegerv(GL_MINOR_VERSION)}; {glGetString(GL_VERSION).decode()}; {glGetString(GL_VENDOR).decode()} -- GL/SL:{glGetString(GL_SHADING_LANGUAGE_VERSION).decode() } -- {glGetString(GL_RENDERER).decode()}\"\n"
pygame.quit()
# Write some general info about the benchmark run.
res_out.write(glnfo)
from datetime import date
today = date.today()
res_out.write(f"\"{today.strftime('%B %d, %Y')}\"\n")
with open(DATA_DIR / "release.txt") as fin:
for line in fin.readlines():
res_out.write(f"\"WolfGPU {line.strip()}\"\n")
res_out.write(f'"Command: {" ".join(sys.argv[:])}"\n')
import platform
if platform.system() == "Windows":
res_out.write(f"\"On {platform.system()}; computer: {os.getenv('COMPUTERNAME')}\"")
else:
res_out.write(f"\"On {platform.system()}\"")
res_out.write("\n")
res_out.write("\"Total Active Tiles: Total number of active tiles, cumulated over the whole simulation.\"\n")
res_out.write("\"Est. T.A.T. : Estimated Total Active Tiles : an estimation of what should be the actual number of active tiles.\"\n")
res_out.write("\"Real/Est. T.A.T. : Ratio of actual over estimated T.A.T.\"\n")
res_out.write("\n")
res_out.write("Nr,Start time,End time,Size,Nb iter,Iterations per second,Active Tiles,Total Time(s),Tiles/s,Total Active Tiles,Est. T.A.T.,Real/Est T.A.T.\n")
res_out.flush()
logging.info(f"Benchmark results will be stored in {BENCHMARK_AGGREGATED_RESULT_FILE.absolute()}")
for sim_ndx, size in enumerate(SIM_SIZES):
# FIXME For the moment, deleting GPU resources doesn't work
# So I force a GPU flush by restarting pygame on each simulation.
pygame.init()
nfo = pygame.display.Info()
window_width, window_height = estimate_best_window_size(nfo.current_w, nfo.current_h, min(SIM_SIZES), min(SIM_SIZES))
pygame.display.set_mode((window_width, window_height), pygame.OPENGL|pygame.DOUBLEBUF, vsync = 0)
pygame.display.set_caption("WolfGPU - Benchmark")
gameIcon = pygame.image.load(
str(DATA_DIR / "wolf_logo.bmp")
)
pygame.display.set_icon(gameIcon)
# Why do we iconify ? Because we have noticed that the performance on
# medium sized problems (1024,2048) are greatly improved when running
# the benchmark on GPU1 over Remote Desktop Connection.
pygame.display.iconify()
print("-"*80)
print(f"Benchmark {sim_ndx+1}/{len(SIM_SIZES)}")
print("-"*80)
cli_args.size = size
# May seem a lot but necessary to have the water to cover the whole simulation surface
# in bigger simulation
cli_args.npas = max(2000,size*3)
m = scenario4(destination_dir, N=int(cli_args.size), active_surface=float(cli_args.ap), optim_pas=True)
glsim = load_sim_to_gpu(m, tile_size=cli_args.tile_size, interpret_freq_type_as_wolf=True, init_gl=False)
force_params_in_gpu(cli_args, glsim)
simulation_runner = SimulationRunner(
glsim,
record_path=destination_dir / "gpu_results",
early_out_delta_max=None,
)
glsim._load_shaders()
now = datetime.now()
start_time = now.strftime("%H:%M:%S")
simulation_runner.full_run(refresh_view=0.5)
now = datetime.now()
end_time = now.strftime("%H:%M:%S")
gs = glsim.read_global_state()
iter_from_time = simulation_runner.iter_per_second * simulation_runner.total_duration
est_active_tiles = float((size**2) / (glsim.tile_size**2)) * iter_from_time
tiles_per_sec = round(gs.total_active_tiles / simulation_runner.total_duration)
res_out.write(
f"{sim_ndx+1},{start_time},{end_time},{size},{cli_args.npas},{simulation_runner.iter_per_second:.2f},{gs.nb_active_tiles},{simulation_runner.total_duration:.1f},{tiles_per_sec},{gs.total_active_tiles},{est_active_tiles:.1f},{gs.total_active_tiles/est_active_tiles:.2f}\n"
)
# Force write to disk. You need the two instructions (see python doc
# for os.fsync())
res_out.flush()
os.fsync(res_out.fileno())
glsim.drop_gl_resources()
from wolfgpu.gl_utils import gl_clear_all_caches
gl_clear_all_caches()
pygame.quit()
logging.info(f"Benchmark file is: {BENCHMARK_AGGREGATED_RESULT_FILE.absolute()}")
[docs]
def main_func():
#from datetime import datetime
global datetime
#datetime.now()
print(LOGOS[datetime.now().second % len(LOGOS)] + "\n")
if (DATA_DIR / "release.txt").exists():
with open(DATA_DIR / "release.txt") as fin:
lines = fin.readlines()
if lines:
print(lines[0].strip())
print("Copyright (c) 2023,2024 Hydraulics in Environmental and ")
print(" Civil Engineering (HECE), University of Liège, Belgium")
if datetime.today().month == 1 and datetime.today().day <= 14:
print("Happy new year !")
print("")
BUILT_IN_SCENARIOS = [
"scenario_cube_drop",
"scenario4",
"small_movement",
"scenario_sine_ground_sine_water",
"scenario_still_water_rocky_bed",
"still_water",
"scenario_multiple_infiltration",
] # , "scenario1", "scenario_drying", "y_steady_influx"]
max_term_width = min(80, shutil.get_terminal_size().columns)
os.environ['COLUMNS'] = str(max_term_width) #shutil.get_terminal_size().columns)
parser = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog= "\n".join(wrap( "When running over a Wolf modelisation, the model is first "
+ "copied to a work directory. That's where you'll get your "
+ "results. We do that because we may run WolfCLI.exe and "
+ "thus override any previous results.", width=max_term_width))
+ "\n\n"
+ "\n".join(wrap( "A bathymetry file is scanned regularly during the simulation. If you "
+ "update it, the simulation will integrate it on its next iteration. The "
+ "file must be a saved numpy array of np.float32, its name must be "
+ "'simul.top_update.npy' and located in the work directory (see info logs).", width=max_term_width)))
parser.add_argument(
"scenario",
nargs="?",
default=None,
help=f"Path to Wolf modelisation directory (we expect 'simul' as generic file) or a built in scenario: {', '.join(BUILT_IN_SCENARIOS)}",
)
parser.add_argument("-o", "-output", default=None, type=str, help="Directory where to store results", metavar="PATH", dest="output")
#parser.add_argument("-help", action="store_true", help="Help")
parser.add_argument(
"-ap",
default=None,
type=float,
help="Wet area proportion (0 to 1) of computation domain (used to parametrize some built in test scenarios)",
)
parser.add_argument(
"-size",
default=None,
type=int,
help="Problem size (in meshes) (used to parametrize some built in test scenarios)",
)
parser.add_argument(
"-skip-write",
action="store_true",
default=False,
help="Skip saving the scenario to disk (only for sy,nthetic scenarios)",
)
parser.add_argument(
"-wolf",
default=None,
type=int,
help="Set number of steps in WolfCLI and start WolfCLI.",
)
parser.add_argument("-npas", default=None, type=int, help="Number of steps to do.")
parser.add_argument("-sim-duration", default=None, type=str,
help="Duration of the simulation (in simulation time). Mutually exclusive with -npas. " \
"The duration is given in seconds by default. It can also be given with a string " \
"such as 1h1m1s for 3661 seconds.")
parser.add_argument(
"-optim-pas",
default=None,
type=int,
help="Force time step optimisation. 0 = no optimisation, 1 = optimisation.",
)
parser.add_argument(
"-dur",
default=None,
type=float,
help="Time step duration (when not optimized).",
)
parser.add_argument(
"-ponderation",
default=None,
help="Force ponderation: 'euler' for single-euler (no ponderation), 0 < PONDERATION < 1 for RungeKutta 2nd order."
" In 2nd order RK, the final equation is PONDERATION*estimator + (1-PONDERATION)*corrector.",
)
parser.add_argument(
"-infilerp",
default=None,
help="Type of infiltration interpolation. Can only be 'linear'. ",
)
parser.add_argument(
"-courant", default=None, type=float, help="Force Courant number."
)
parser.add_argument(
"-froude-max", default=None, type=float, help="Force Froude maximum."
)
parser.add_argument(
"-manning", default=None, type=float, help="Force Manning coefficient."
)
parser.add_argument(
"-freq", type=str, default=None, help="If an integer number, one record every 'freq' iteration. If an integer number N suffixed by 's' (e.g. '10s'), one record every N seconds (measured in simulation time). "
)
parser.add_argument(
"-tile-size", default=16, type=int, help="tile size, must be a power of 2"
)
parser.add_argument(
"-dry-up-loops",
type=int,
default=None, # Default value used when an argument is not provided
metavar="N",
help="Exact number of dry up outer loops (0 means no dry up management, 1 means one pass, etc.). Outer loops are those used to solve dry ups across tiles. If you don't provide this parameter, then we will do as many loops as necessary to solve all the dry ups.",
)
# parser.add_argument(
# "-dry-up",
# default="exact",
# type=str,
# help="Select the dry up management strategy",
# )
parser.add_argument(
"-no-tiles-packing",
default=False,
action="store_true",
help="Disable tiles packing",
)
parser.add_argument(
"-early-out-threshold",
default=None,
type=float,
help="If specified, will stop simulation if for 3 records in a row, sum{all meshes m}(abs(m.h) + abs(m.qx) + abs(m.qy)) <= EARLY_OUT_THRESHOLD",
)
parser.add_argument(
"-start-from-last-record",
default=False,
action="store_true",
help="Restart the sim from the last record",
)
parser.add_argument(
"-restart-from-record",
type=int,
help="Restart the sim from the n-th record (we insist: a record, not a step)",
)
parser.add_argument(
"-threshold-depth",
nargs="?",
const=1e-5,
type=float,
help="Clear all meshes where water depth is below the given level before starting simulation (default 1e-5m)",
)
parser.add_argument(
"-timeline",
nargs=2,
type=str,
help="Plot the evolution of the simulation as reported in the results. First parameter is the path to the work directory (containing the Wolf model), second parameter is the path to the results (both as reported when starting the simulation).",
)
parser.add_argument(
"-benchmark",
default=False,
action="store_true",
help=f"Run a full benchmark and report the results in a file called \"{BENCHMARK_AGGREGATED_RESULT_FILE}\".",
)
parser.add_argument(
"-debug", action="store_true", default=False, help="Activate debug logging"
)
parser.add_argument(
"-record-alphas",
action="store_true",
default=False,
help="Enable the recording of the dry up management texture (a.k.a. alpha textures)",
)
parser.add_argument(
"-version",
action="store_true",
default=False,
help="Show information about this program",
)
parser.add_argument(
"-gpu-info",
action="store_true",
default=False,
help="Show information about the GPU",
)
parser.add_argument(
"-quickrun",
type=str,
default="",
metavar="PATH",
help="Run a quick simulation from a SimpleSimulation object",
)
parser.add_argument(
"-no-validation",
action="store_true",
default=False,
help="If set, failed checks made to ensure a simulation is valid will produce warnings instead of errors.",
)
parser.add_argument(
"-optimize-indirection",
action="store_true",
default=False,
help="Activate tile indirection optimization for faster wet tiles detection (still being tested).",
)
if POWER_MODE:
parser.add_argument(
"-shader-log",
type=Path,
default=None,
metavar="PATH",
help="(PowerMode) Enable logging of shader codegen to the given path.",
)
parser.add_argument(
"-profiler",
default=False,
action="store_true",
help="Activate profiling",
)
# parser.add_argument('-optim-tiles', choices=["CS","RS"], default=None, help="Enable or disable tile geometry shader GPU optimisation. CS=Compute shader, RS=Regular shader.")
# Transform double hyphen in simple hyphen (allows one to use
# Unix like arguments.
for i in range(len(sys.argv)):
if sys.argv[i] != "--help" and sys.argv[i].startswith("--"):
sys.argv[i] = sys.argv[i][1:]
cli_args = parser.parse_args()
if cli_args.debug:
init_global_logging(logging.DEBUG)
else:
init_global_logging(logging.INFO)
assert not cli_args.wolf or not bool(
cli_args.start_from_last_record
), "-npas-wolf incompatible with -start-from-last-record"
assert not cli_args.wolf or not bool(
cli_args.restart_from_record
), "-npas-wolf incompatible with -restart-from-record"
assert not (
cli_args.start_from_last_record and cli_args.restart_from_record
), f"-start-from-last-record and -restart-from-record options are mutually exsclusive."
assert (cli_args.sim_duration is None and cli_args.npas is None) or ((cli_args.sim_duration is not None) ^ (cli_args.npas is not None)), \
"-sim-duration and -npas parameters are mutually exclusive."
assert not cli_args.optimize_indirection or (cli_args.optimize_indirection and not cli_args.no_tiles_packing), \
"For optimized indirection, you must have tile packing"
if cli_args.scenario and cli_args.quickrun:
logging.error(
f"Using -quick-run and passing a scenario ({cli_args.scenario}) on the command line doesn't make sense. Choose one or the other."
)
exit()
if cli_args.scenario in BUILT_IN_SCENARIOS:
if cli_args.sim_duration is not None and "s" in cli_args.sim_duration.lower():
logging.warning("The -sim-duration parameters set with time will not be written in the Wolf simulation data.")
if cli_args.ponderation is not None:
try:
cli_args.ponderation = float(cli_args.ponderation)
assert (
0.0 < float(cli_args.ponderation) < 1.0
), "Ponderation value is not correct. It should be between 0 and 1"
except:
# Not a float, so a string.
cli_args.ponderation = cli_args.ponderation.lower()
assert (
cli_args.ponderation == "euler"
), f"I don't recognize the ponderation : {cli_args.ponderation}"
wolf_sim_path = cli_args.ponderation
if cli_args.gpu_info:
from gl_utils import query_gl_caps, init_gl
init_gl(10, 10)
query_gl_caps()
exit()
if cli_args.version:
with open(DATA_DIR / "release.txt") as fin:
print(fin.readlines()[1])
exit()
if cli_args.benchmark:
benchmark(Path.home() / "test_sim", cli_args)
exit()
if cli_args.timeline:
plot_evolution_of_simulation(
Path(cli_args.timeline[0]), Path(cli_args.timeline[1])
)
exit()
if POWER_MODE and cli_args.scenario is None and not cli_args.quickrun:
# FIXME Dirty shortcuts to develop faster (should put that in VSCode
# .json config, but it's painful).
cli_args.scenario = r"d:\StephaneC\Theux"
if cli_args.scenario is None and not cli_args.quickrun:
parser.print_help()
exit()
simulation_dir = None
if cli_args.scenario in BUILT_IN_SCENARIOS:
simulation_dir = (Path.home() / "test_sim" / cli_args.scenario).resolve()
elif cli_args.scenario:
simulation_dir = Path(cli_args.scenario).resolve()
elif cli_args.quickrun:
simulation_dir = Path(cli_args.quickrun).resolve()
else:
logging.error("Could not figure out simulation directory")
exit()
if cli_args.output:
destination_dir = Path(cli_args.output).resolve()
if destination_dir.exists():
assert not destination_dir.is_file(), f"The output directory ({destination_dir}) is a file ???"
assert destination_dir.is_dir(), f"The output directory ({destination_dir}) is not a directory ?"
assert simulation_dir.resolve() != destination_dir.resolve(), f"The result directory {simulation_dir} is the same as the model directory."
else:
destination_dir = simulation_dir / "simul_gpu_results"
glsim = None
if cli_args.scenario in BUILT_IN_SCENARIOS:
logging.info("Building test scenario")
if cli_args.ponderation is None:
cli_args.ponderation = 0.3
if cli_args.courant is None:
cli_args.courant = 0.25
if cli_args.manning is None:
cli_args.manning = 0.04
if cli_args.scenario == "still_water":
assert (
cli_args.ap is None
), "-ap is of no use for the cube drop scenario. Only -size is."
if cli_args.size is None:
cli_args.size = 128
# Perfectly still water
m = scenario_still_water(
simulation_dir,
width=cli_args.size,
height=cli_args.size,
rockyness=0,
water_height=np.float32(5.0),
)
elif cli_args.scenario == "scenario_cube_drop":
assert (
cli_args.ap is None
), "-ap is of no use for the cube drop scenario. Only -size is."
if cli_args.size is None:
cli_args.size = 128
m = scenario_cube_drop(
simulation_dir, width=cli_args.size, height=cli_args.size
)
elif cli_args.scenario == "scenario_sine_ground_sine_water":
assert (
cli_args.ap is None
), "-ap parameter is not used for scenario_sine_ground_sine_water"
if cli_args.size is not None:
assert cli_args.size >= 100, "-size parameter must be >= 100"
else:
cli_args.size = 100
m = scenario_sine_ground_sine_water(simulation_dir, N=cli_args.size)
elif cli_args.scenario == "scenario_still_water_rocky_bed":
assert (
cli_args.ap is None
), "-ap parameter is not used for scenario_sine_ground_sine_water"
if cli_args.size is not None:
assert cli_args.size >= 10, "-size parameter must be >= 10"
else:
cli_args.size = 10
m = scenario_still_water_rocky_bed(
simulation_dir, width=cli_args.size, height=cli_args.size
)
elif cli_args.scenario == "scenario4":
if cli_args.ap is None:
cli_args.ap = 0.8
if cli_args.size is None:
cli_args.size = 128
assert (
0.5 <= cli_args.ap < 1.0
), "Proportion of occupied meshes (-ap x) must me between >= 0.5 and < 1.0"
assert (
cli_args.size >= 16
), f"Size of problem (param. -size) must be >= 16, it is {cli_args.size}"
m = scenario4(
simulation_dir,
N=int(cli_args.size),
active_surface=float(cli_args.ap),
optim_pas=True,
)
elif cli_args.scenario == "small_movement":
assert cli_args.ap is None, "-ap parameter is useless in this scenario"
assert cli_args.size is None, "-size parameter is useless in this scenario"
m = scenario_small_movement(simulation_dir)
m.parameters._nb_timesteps = 300_000 # OK for rust, not on my PC
m.parameters._scheme_optimize_timestep = 1
m.parameters._scheme_rk = 0.3
# Force some parameter
cli_args.ponderation = m.parameters._scheme_rk
cli_args.courant = m.parameters._scheme_cfl
cli_args.manning = 0.04
elif cli_args.scenario == "scenario1":
raise Exception(
"This scenario uses boundary conditions which are still disabled"
)
assert cli_args.size >= 32, "Size of problem (-size) must be >= 32"
m = scenario1(
simulation_dir,
size=int(cli_args.size),
bc_qx_factor=1.0,
bc_hmod_factor=1.0,
WATER_HEIGHT=10.0,
IC_SPEED=1.0,
friction=0.0,
)
force_cli_params_in_wolf_model(cli_args, m)
write_simulation(m, WOLF_CLI_EXE)
elif cli_args.scenario == "scenario_drying":
m = scenario_drying(
simulation_dir,
size=int(cli_args.size),
WATER_HEIGHT=0.2,
friction=0.0,
top_topo=1.0,
)
# import matplotlib.pyplot as plt
# plt.imshow(m.top.array[2:-2,2:-2])
# plt.show()
m.parameters._nb_timesteps = int(cli_args.npas)
m.parameters._scheme_optimize_timestep = int(cli_args.optim_pas)
m.parameters._timestep_duration = float(cli_args.dur)
m.parameters._writing_frequency = cli_args.freq
m.parameters._scheme_k_venkatakrishnan = 0.0 # was 1.0
m.parameters._scheme_centered_slope = 1 # was 2
write_simulation(m, WOLF_CLI_EXE)
elif cli_args.scenario == "y_steady_influx":
m = create_steady_flow_swimming_pool_bottom_to_top(
simulation_dir, width=cli_args.size // 2, height=cli_args.size
)
elif cli_args.scenario in BUILT_IN_SCENARIOS:
m = locals()[cli_args.scenario](simulation_dir)
elif not cli_args.quickrun:
# User gives us a path to a regular, existing Wolf simulation
wolf_sim_path = Path(cli_args.scenario)
if not wolf_sim_path.exists():
logging.error(f"I can't find the scenario '{cli_args.scenario}'. I have looked in: {wolf_sim_path.absolute()}")
exit()
wolf_sim_path = wolf_sim_path.resolve()
# Loading a regular Wolf modelisation
assert (
cli_args.ap is None
), f"Setting -ap doesn't make sense when loading a Wolf scenario"
assert (
cli_args.size is None
), f"Setting -size doesn't make sense when loading a Wolf scenario"
if wolf_sim_path.is_file():
generic_file = wolf_sim_path.name
wolf_sim_path = wolf_sim_path.parent
else:
generic_file = "simul"
assert (
wolf_sim_path / generic_file
).exists(), f"I can't find {wolf_sim_path / generic_file}"
logging.warning(
"Assuming 'simul' is the name of the generic file in the Wolf directory"
)
# When we use a built in scenario we write it as a WolfModel first
# to be able to use the wolf model loader afterewards (to test it)
if cli_args.scenario in BUILT_IN_SCENARIOS:
if not cli_args.skip_write:
write_simulation(m, WOLF_CLI_EXE)
generic_file = "simul"
# Load the wold modelisation in the GPU
logging.info(
f"Loading synthetic model at {str(simulation_dir / generic_file)}"
)
m = prev_sim2D( str(simulation_dir / generic_file))
elif cli_args.scenario:
logging.info(f"Loading Wolf model at {str(Path(cli_args.scenario) / generic_file)}")
m = prev_sim2D( str(Path(cli_args.scenario) / generic_file) )
if cli_args.scenario:
force_cli_params_in_wolf_model(cli_args, m)
m.force_load()
if cli_args.threshold_depth:
candidates = m.hbin.array[:, :] > 0
msk = (m.hbin.array[:, :] < cli_args.threshold_depth) & candidates
# from matplotlib import pyplot as plt
# m.hbin.array[msk] = 1000
# plt.imshow(m.hbin.array)
# plt.show()
m.hbin.array[msk] = 0
logging.info(
f"Thresholding water depth < {cli_args.threshold_depth} m. {np.sum(msk)} ({np.sum(msk)/np.sum(candidates):.3f}%) meshes were zeroed."
)
if POWER_MODE and (cli_args.shader_log is not None and not cli_args.shader_log.exists()):
cli_args.shader_log.mkdir(parents=True)
shader_log_path = cli_args.shader_log
else:
shader_log_path = None
if cli_args.no_tiles_packing:
tiles_packing_mode = TilePackingMode.TRANSPARENT
else:
tiles_packing_mode = TilePackingMode.REGULAR
if cli_args.scenario in BUILT_IN_SCENARIOS:
if cli_args.scenario == "scenario4":
report_label = f"scenario4_portion_{int(float(cli_args.ap)*10)}"
else:
report_label = re.sub(
"_+",
"_",
Path(cli_args.scenario).name.replace(" ", "_").replace("-", "_"),
)
else:
report_label = "report_sim"
if cli_args.quickrun:
if not Path(cli_args.quickrun).exists():
logging.error("You must give an existing SimpleSim directory to -quickrun")
exit()
sim_path = Path(cli_args.quickrun).resolve()
logging.info(f"Loading simple simulation {sim_path}")
simple_simulation = SimpleSimulation.load(sim_path)
if cli_args.manning is not None:
simple_simulation.manning[:, :] = cli_args.manning
page_flipper = init_pygame(simple_simulation.param_nx, simple_simulation.param_ny)
glsim = load_simple_sim_to_gpu(
simple_simulation,
tile_size=cli_args.tile_size,
shader_log_path=shader_log_path,
tiles_packing_mode=tiles_packing_mode,
optimize_indirection=cli_args.optimize_indirection,
fail_if_invalid_sim=not cli_args.no_validation
)
else:
page_flipper = init_pygame(m.parameters.nbx, m.parameters.nby)
glsim = load_sim_to_gpu(
m, tile_size=cli_args.tile_size, shader_log_path=shader_log_path, interpret_freq_type_as_wolf=True,
tiles_packing_mode=tiles_packing_mode,
optimize_indirection=cli_args.optimize_indirection,
init_gl=False,
fail_if_invalid_sim=not cli_args.no_validation
)
glsim.set_optim_geom_shaders(TileOptimisation.COMPUTE_SHADER)
force_params_in_gpu(cli_args, glsim)
pygame.display.set_caption("WolfGPU")
gameIcon = pygame.image.load(
str(DATA_DIR / "wolf_logo.bmp")
)
pygame.display.set_icon(gameIcon)
# if not cli_args.skip_write:
# if cli_args.scenario in BUILT_IN_SCENARIOS:
# short_normat_name = cli_args.scenario
# elif cli_args.scenario:
# short_normat_name = wolf_sim_path.name
# save_short_format(m, glsim, simulation_dir / short_normat_name)
# FIXME Rename to GPUSimulationRunner
simulation_runner = SimulationRunner(
glsim,
record_path=destination_dir,
early_out_delta_max=cli_args.early_out_threshold,
page_flip_func=page_flipper,
enable_alpha_recording=cli_args.record_alphas
)
if not cli_args.start_from_last_record and not cli_args.restart_from_record:
res_dir = Path(simulation_runner.record_file())
logging.info(f"Clearing the results directory {res_dir}")
if os.path.exists(res_dir):
delete_dir_recursion(res_dir)
# shutil.rmtree(destination_dir)
if not os.path.exists(res_dir):
os.makedirs(res_dir)
# from tests.arrayview import array_view
# array_view(m.napbin.array.data)
if cli_args.quickrun:
print_sim_params(glsim, simulation_runner, simple_simulation.manning, simple_simulation.nap, simple_simulation.h)
else:
print_sim_params(glsim, simulation_runner, m.frot.array.data, m.napbin.array.data, m.hbin.array.data)
logging.info(f"GPU Results directory : {simulation_runner.record_file()}")
logging.info(f"Bathymetry update file is scanned at {simulation_runner.scanned_file()}")
print()
logging.info(
"Simulating... Hit 'q' on the simulation view to abort (ctrl-c is no good)"
)
# cProfile.run('simulation_runner.full_run(refresh_view=0.5)', 'stats')
if cli_args.start_from_last_record:
simulation_runner.restart_from_record()
elif cli_args.restart_from_record:
simulation_runner.restart_from_record(cli_args.restart_from_record)
if (
cli_args.start_from_last_record or cli_args.restart_from_record
) and cli_args.npas:
assert cli_args.npas > simulation_runner.step_num, (
"The -npas must be the *total* number of steps in the "
"simulation (not a number of steps you want to add on top of what's "
"already done). The simulation you pointed at currently has "
f"{simulation_runner.step_num} iterations."
)
if POWER_MODE and cli_args.profiler:
logging.info("Profiler is on")
with cProfile.Profile() as pr:
simulation_runner.full_run(refresh_view=0.1, label=f"sim_{report_label}.pdf")
pr.dump_stats("stats")
else:
logging.info("Profiler is off")
simulation_runner.full_run(refresh_view=0.1, label=f"sim_{report_label}.pdf")
# from gl_utils import read_texture, GL_COLOR_ATTACHMENT0, GL_R32I
# print(read_texture(simulation_runner._glsim._infiltration_fb, GL_COLOR_ATTACHMENT0, simulation_runner._glsim.width, simulation_runner._glsim.height, GL_R32I))
pygame.display.quit()
if cli_args.wolf is not None:
# Clean tqdm's output
# Maybe useless: https://github.com/tqdm/tqdm/issues/737
# sys.stdout.flush()
# print("")
# Prepare the wolf modelisation to be run by WolfCLI
generic_path = Path(m.mydir) / Path(m.filenamegen)
m.parameters._nb_timesteps = cli_args.wolf
m.parameters.write_file(fn=str(generic_path))
# # FIXME VERY BIG HACK See
# # https://gitlab.uliege.be/HECE/HECEPython/-/issues/38
# with open(str(generic_path) + ".par", "a") as f:
# for b in range(1, m.myparam.nblocks + 1):
# f.write(f"{b}\n")
# for b in range(m.myparam.nblocks):
# f.write(f'""\n')
wolf_cli_process = run_wolf_cli(Path(m.mydir))
REG = re.compile(r"\*s\*3 ITERATION : +([0-9]+)")
REG_CPU = re.compile(r"\*s\*1 TEMPS TOTAL CPU : +([0-9.]+(E-[0-9]+)?)")
last_wolfcli_iter = None
logging.info(f"Watching {Path(destination_dir)/'simul.NFO'}")
# model: prev_sim2D = prev_sim2D(dir=str(destination_dir), splash=False, in_gui=False)
# sys.stdout trick, see : https://stackoverflow.com/questions/45862715/how-to-flush-tqdm-progress-bar-explicitly
progress_bar = tqdm.tqdm(range(cli_args.wolf), leave=False, file=sys.stdout)
nb_errors = 0
wolf_nfo = Path(destination_dir) / "simul.NFO"
while wolf_cli_process.poll() is None:
try:
if wolf_nfo.exists() and wolf_nfo.stat().st_size > 1024:
with open(Path(destination_dir) / "simul.NFO", "rb") as file:
file.seek(-1024, os.SEEK_END)
s = file.read().decode()
match = REG.findall(s)
if match:
i = match[-1]
progress_bar.update(int(i) - progress_bar.n)
except Exception as ex:
logging.error(ex)
nb_errors += 1
if nb_errors >= 5:
logging.error(
"Too many errors while watching WolfCLI.exe running. Did wolfcli.exe accept the modelisation ? Aborting."
)
wolf_cli_process.kill()
exit()
sleep(1)
# Clean tqdm's output
sys.stdout.flush()
print("")
with open(Path(destination_dir) / "simul.NFO", "rb") as file:
nfo = file.read().decode(encoding="ISO-8859-15")
match = REG_CPU.findall(nfo)
assert match, "I'm expecting something in wolfcli.exe's .NFO file"
# Remove the setup time of WolfCLI
cpu_time = float(match[-1][0]) - float(match[0][0])
with open("report.txt", "a") as fout:
gpu_iter_per_sec = max(1, simulation_runner.iter_per_second)
wolfcli_iter_per_sec = max(1, float(cli_args.wolf) / cpu_time)
accel = gpu_iter_per_sec / wolfcli_iter_per_sec
wolf_sim_path = f"Total time CPU WolfCLI {cpu_time:.3f}s for {cli_args.wolf} iterations. => {wolfcli_iter_per_sec:g} it/s. GPU = {gpu_iter_per_sec:.2f} it/s. Accel={accel:.1f}\n"
fout.write(wolf_sim_path)
logging.info(wolf_sim_path)
""" FIXME Now the question is : how do I run that
in a way that I don't block WolfGUI ?
- threads: the GIL will block anything CPU bound
and I wonder if OpenGL counts as IO bound...
- multiprocess: fine but howdir do I communictate
progress information ? can I access the OpenGL
context ?
"""
if __name__ == "__main__":
main_func()