diff --git a/.coveragerc b/.coveragerc index 9adae452..a7cba45f 100755 --- a/.coveragerc +++ b/.coveragerc @@ -1,4 +1,3 @@ [run] omit = - src/lisflood/hydrological_modules/compile_kinematic_wave_parallel_tools.py tests/* diff --git a/.github/workflows/ci_env.yml b/.github/workflows/ci_env.yml new file mode 100644 index 00000000..0e2944bd --- /dev/null +++ b/.github/workflows/ci_env.yml @@ -0,0 +1,47 @@ +name: Lisflood OS Unit Tests + +on: [push] + +jobs: + tests: + runs-on: ubuntu-20.04 + strategy: + max-parallel: 5 + + steps: + - uses: actions/checkout@v3 + - uses: conda-incubator/setup-miniconda@v2 + with: + auto-update-conda: true + python-version: 3.7 + - name: Conda info + shell: bash -el {0} + run: conda info + - name: Install python and gcc + shell: bash -el {0} + run: | + conda install -c conda-forge python=3.7 + conda install -c conda-forge gcc=12.1.0 + - name: Install gdal and pcraster + shell: bash -el {0} + run: | + conda install -c conda-forge gdal pcraster + - name: Install dependencies + shell: bash -el {0} + run: | + pip install -r requirements.txt + - name: Install lisflood-module + shell: bash -el {0} + run: | + pip install . + - name: Check installation + shell: bash -el {0} + run: | + gdal-config --version + python -c "from osgeo import gdal; print(gdal.__version__)" + conda list + - name: Test with pytest + shell: bash -el {0} + run: | + pip install pytest pytest-cov + pytest \ No newline at end of file diff --git a/.gitignore b/.gitignore index 68d80100..222e9834 100644 --- a/.gitignore +++ b/.gitignore @@ -17,13 +17,11 @@ npm-debug.log* __pycache__/ build/ *.so -kinematic_wave_parallel_tools.html lisflood_model.egg-info /dist/ .eggs .tox /src/dist/ -/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.c .coverage *.tox.ini .vscode/ diff --git a/Dockerfile b/Dockerfile index b6c2bb7d..92ad8fb4 100755 --- a/Dockerfile +++ b/Dockerfile @@ -27,9 +27,6 @@ COPY src/lisfloodSettings_reference.xml / COPY LICENSE / COPY VERSION / -# Compile kwpt -RUN cd /lisflood/hydrological_modules && conda run -n lisflood python compile_kinematic_wave_parallel_tools.py build_ext --inplace - # RUN Tests COPY tests/. /tests/ COPY pytest.ini /tests diff --git a/MANIFEST.in b/MANIFEST.in index f83ba9a5..bb1d0bd5 100755 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,4 +1,3 @@ -include src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.pyx include requirements.txt include *.xml include VERSION diff --git a/README.md b/README.md index 79582933..8441f971 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,7 @@ You can use conda environment to easily install dependencies. ```bash conda create --name lisflood python=3.7 -c conda-forge conda activate lisflood -conda install -c conda-forge pcraster +conda install -c conda-forge pcraster gdal ``` * Install lisflood-model pypi package @@ -61,26 +61,7 @@ pip install -r requirements.txt If you don't use conda but a plain virtualenv, you need to install PCRaster and GDAL by your own and include its python interface in PYTHONPATH environment variable. For details, please follow instruction on [official docs](http://pcraster.geo.uu.nl). - -3. Compile the cython module kinematic_wave_parallel_tool - -To compile this Cython module to enable OpenMP multithreading (parallel kinematic wave): - -* Delete the files *.so (if any) in directory hydrological-modules - -* Inside the hydrological_modules folder, execute "python compile_kinematic_wave_parallel_tools.py build_ext --inplace" - -Important: the module has to be compiled on the machine where the model is run - the resulting binary is not portable. - -Then in the settings file the option "numberParallelThreadsKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> -```xml - -``` -4. Run a cold run for the test catchment +3. Run a cold run for the test catchment Now your environment should be set up to run lisflood. Try with a prepared settings file for one of the two test catchments: diff --git a/VERSION b/VERSION index ef8d7569..81911389 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -4.2.0 \ No newline at end of file +4.3.0 \ No newline at end of file diff --git a/docs/3_step2_installation/index.md b/docs/3_step2_installation/index.md index ced4871a..db2054ec 100644 --- a/docs/3_step2_installation/index.md +++ b/docs/3_step2_installation/index.md @@ -95,32 +95,12 @@ If you don't use conda but a plain virtualenv, you need to install PCRaster and For details, please follow instruction on [official docs](https://pcraster.geo.uu.nl/pcraster/4.3.1/documentation/pcraster_project/install.html). -3. **Compile the cython module kinematic_wave_parallel_tool** - - To compile this Cython module to enable OpenMP multithreading (parallel kinematic wave): - - * Delete the files *.so (if any) in directory hydrological-modules - * Inside the hydrological_modules folder, execute `python compile_kinematic_wave_parallel_tools.py build_ext --inplace` - - Important: the module has to be compiled on the machine where the model is run - the resulting binary is not portable. - - Then in the settings file the option "numberParallelThreadsKinematicWave" may take the following values: - - * "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) (do not set it if other simulations are running on the same machine/node) - * "1" : serial execution (not parallel) - * "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") - - ```xml - - ``` - -4. **Run a cold run for the test catchment** +3. **Run a cold run for the test catchment** Now that your environment should be set up to run lisflood, you may try with a prepared settings file for test catchment included into the tests/data folder: ```bash - python src/lisf1.py tests/data//settings/cold_day_base.xml + python src/lisf1.py tests/data//settings/cold.xml ``` 4. **Run LISFLOOD unit tests** diff --git a/docs/3_step3_preparing-setting-file/lisfloodSettings_reference.xml b/docs/3_step3_preparing-setting-file/lisfloodSettings_reference.xml index 8021c269..072b9a98 100644 --- a/docs/3_step3_preparing-setting-file/lisfloodSettings_reference.xml +++ b/docs/3_step3_preparing-setting-file/lisfloodSettings_reference.xml @@ -40,6 +40,7 @@ You can use builtin path variables in this template and reference to other paths + @@ -93,6 +94,9 @@ You can use builtin path variables in this template and reference to other paths # report maps + + + @@ -102,10 +106,10 @@ You can use builtin path variables in this template and reference to other paths + - - - + + @@ -117,20 +121,48 @@ You can use builtin path variables in this template and reference to other paths +************************************************************** +netCDF parameters +************************************************************** + +!-- Optimization of netCDF I/O through chunking and caching. + +The option "NetCDFTimeChunks" may take the following values: + - "-1" : Load everything upfront + - "[positive integer number]" : Chunk size defined by user + - "auto" : Let xarray to decide + + + + + +The option "MapsCaching" may take the following values: + - "True" : Cache maps during execution + - "False" : No caching + + + + + +The option "OutputMapsChunks" may take the following values: + - "[positive integer number]" : Dump outputs to disk every X steps (default 1) + + + ************************************************************** -PARALLELIZED KINEMATIC ROUTING +PARALLELISATION WITH NUMBA (USED IN ROUTING AND SOILLOOP) ************************************************************** -!-- Parallelisation of the kinematic wave routing (via openMP multi-threading). -The option "numberParallelThreadsKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) - (do not set it if other simulations are running on the same machine/node) +!-- Parallelisation of using Numba runtime compiled library . +The option "numCPUs_parallelNumba" may take the following values: + - "0" : set to NUMBA_NUM_THREADS Environment Variable + (if NUMBA_NUM_THREADS is not set, will take the number of CPU cores determined by python's multiprocessing.cpu_count()) - "1" : serial execution (not parallel) - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> - + (if exceeding NUMBA_NUM_THREADS, the value is set to NUMBA_NUM_THREADS) --> + ************************************************************** @@ -490,7 +522,7 @@ Path of the initial value maps e.g. lzavin.map (org=$(PathRoot)/outPo) - + Static maps path /perm/mo/mocm/proj/efas/xdom/data/staticData/lisflood @@ -1105,6 +1137,18 @@ More than the transpiration is added e.g to prevent salinisation + + +Annual amount (m3) of treated wastewater reused for irrigation + + + + + +Number of days over which the annual amount of treated wastewater for irrigation is used + + + Consumptive Use (1-Recycling ratio) for livestock water use (0-1) [-] @@ -1573,6 +1617,65 @@ otherwise doesn't influence model results at all) + +PF MAPS + + + + +Reported pF soil layer 1a, other fraction + + + + + +Reported pF soil layer 1b, other fraction + + + + + +Reported pF soil layer 2, other fraction + + + + + +Reported pF soil layer 1a, forest fraction + + + + + +Reported pF soil layer 1b, forest fraction + + + + + +Reported pF layer 2, forest fraction + + + + + +Reported pF soil layer 1a, irrigation fraction + + + + + +Reported pF soil layer 1b, irrigation fraction + + + + + +Reported pF soil layer 2, irrigation fraction + + + + @@ -1591,20 +1694,8 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT ************************************************************** --> - - - - -Chunking in time dimension used when reading NetCDF forcings (or other NetCDF temporal inputs). -Special values: -"auto": recommended, uses default XArray chunking (depends on the size and type of the array) -"-1": load the whole time series in memory (fastest option but memory consuming) - + - - -Option designed for the calibration. All the static maps and forcings will be stored in a cache so that they don't have to be loaded by each lisflood instance. This option sets the value of NetCDFTimeChunks to -1, meaning that the whole time series in the NetCDF inputs is loaded into memory. - @@ -1612,6 +1703,10 @@ This is necessary when using projected coordinates (x,y) + + + + Clone map @@ -2561,6 +2656,18 @@ More than the transpiration is added e.g to prevent salinisation + + +Annual amount (m3) of treated wastewater reused for irrigation + + + + + +Number of days over which the annual amount of treated wastewater for irrigation is used + + + Irrigation crop coefficient @@ -2930,6 +3037,30 @@ Reported storage in upper response box [mm] + + +Reported transpiration maps, other fraction [mm] + + + + + +Reported transpiration maps, forest fraction [mm] + + + + + +Reported transpiration maps, irrigated fraction [mm] + + + + + +Reported transpiration maps, weighted sum over the fractions of each pixel [mm] + + + @@ -3459,15 +3590,63 @@ Reported preferential flow [mm] - + + +Reported percolation from soil layer 1a to soil layer 1b, other fraction [mm] + + + + + +Reported percolation from soil layer 1a to soil layer 1b, forest fraction [mm] + + + + + +Reported percolation from soil layer 1a to soil layer 1b, irrigation fraction [mm] + + + + + +Reported percolation from soil layer 1b to soil layer 2, other fraction [mm] + + + + + +Reported percolation from soil layer 1b to soil layer 2, forest fraction [mm] + + + + + +Reported percolation from soil layer 1b to soil layer 2, irrigation fraction [mm] + + + + + +Reported seepage to groundwater(from soil layer 2 to groundwater), other fraction [mm] + + + + + +Reported seepage to groundwater(from soil layer 2 to groundwater), forest fraction [mm] + + + + -Reported percolation from soil layer 1 to soil layer 2 [mm] +Reported seepage to groundwater(from soil layer 2 to groundwater), irrigation fraction [mm] - + -Reported seepage to groundwater(from soil layer 2 to groundwater) [mm] +Reported seepage to groundwater(from soil layer 2 to groundwater), weighted sum over the fraction of each pixel [mm] @@ -3764,19 +3943,25 @@ Water stored as interception [mm] -Soil moisture upper layer [mm3/mm3] +Soil moisture layer 1a [cu mm / cu mm] -Soil moisture lower layer [mm3/mm3] +Soil moisture lower layer [cu mm/cu mm] -Soil moisture lower layer [mm3/mm3] +Soil moisture layer 1b [cu mm / cu mm] + + + + + +Soil moisture layer 2 [cu mm / cu mm] @@ -4159,6 +4344,30 @@ Reported mass balance error at outlet (as mm water slice) + + +Reported mass balance error due to the split routing module, for each catchment [m3] + + + + + +Reported discharge error at the outlet as a consequence of the mass balance error within the split routing module, for each catchment [m3/s] + + + + + +Reported ratio between the mass balance error and the water volume storage, for each catchment [m3/m3] + + + + + +Average value of the sum of the fractions within a catchment: this value must be 1.0 in order to avoid mass balace errors! [-] + + + Number of cells with Theta lower than @@ -5306,15 +5515,9 @@ TotalAbstractionFromSurfaceWater [mm] summed up for water regions - - -TotalAbstractionFromGroundwater [mm] summed up for water regions - - - - + -TotalIrrigationAbstraction [m3] +Total Abstraction From SurfaceWater and Groundwater [mm] summed up for water regions, montly values @@ -5324,55 +5527,6 @@ TotalPaddyRiceIrrigationAbstraction [m3] - - -TotalLivestockAbstraction [m3] - - - - - -DomesticConsumptiveUseMM [mm] - - - - - -LivestockConsumptiveUseMM [mm] - - - - - -IndustrialConsumptiveUseMM [mm] - - - - - -EnergyConsumptiveUseMM [mm] - - - - - -IrrigationWaterAbstraction [mm] - - - - - - -Water use (water demand is reduced by water available) - - - - - -Water use (water demand is reduced by water available) summed up for water regions - - - Water Exploitation Index - Consumption; for water regions (this is the official EU WEI+): consumption / internal and external availability @@ -5391,6 +5545,11 @@ Water Exploitation Index - Demand; for water regions: demand / internal and exte + + +Water Exploitation Index - Consumption; for water regions + + @@ -5404,21 +5563,15 @@ External available water - + Region Water Consumption - + -Region Water Demand - - - - - -Region Water Demand +Region Water Consumption @@ -5520,12 +5673,6 @@ Abstraction from surface water in m3 per timestep - - -Region abstraction from groundwater in m3 per timestep - - - Region abstraction from surface water in m3 per timestep @@ -5562,24 +5709,6 @@ LakeStorage in M3 - - -Total Abstraction Demand - - - - - -WUseAddM3 - - - - - -totalAddM3 - - - AreatotalIrrigationUseM3 @@ -5598,6 +5727,36 @@ EFlowIndicator (1 on day with ChanQ smaller than EflowThreshold, 0 on normal day + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -5646,4 +5805,3 @@ ReportSteps=$(ReportSteps); - diff --git a/docs/3_step7_commmand-line-flags/index.md b/docs/3_step7_commmand-line-flags/index.md new file mode 100644 index 00000000..8e94474d --- /dev/null +++ b/docs/3_step7_commmand-line-flags/index.md @@ -0,0 +1,43 @@ +# Step 7: Command line flags +LISFLOOD command line takes the following flags as additional arguments after the xml setting file: + +``` + -q --quiet output progression given as . + -v --veryquiet no output progression is given + -l --loud output progression given as time step, date and discharge + -c --checkfiles input maps and stack maps are checked, output for each input map BUT no model run + -n --nancheck check input maps and routing output for any NaN value generated during model run + -h --noheader .tss file have no header and start immediately with the time series + -d --debug debug outputs + -i --initonly only run initialisation, not the dynamic loop + -s --skipvalreplace skip replacement of invalid values in meteo input maps (ignore valid_min and valid_max) +``` + +The flags are utility flags and do not change the behaviour/parameters of the model. Here are the operational details of each flag. + +- **-q --quiet output progression given as .** + The default on-screen output of the lisflood command is the step count and the date/time of each computational step (each step being the run of all the activated modules in sequence for each time step). By setting this "-q" flag only a dot "." will be writtend on stdout for each computational step. + +- **-v --veryquiet no output progression is given** + The default on-screen output of the lisflood command is the step count and the date/time of each computational step (each step being the run of all the activated modules in sequence for each time step). By setting this "-v" flag there will be no output written on stdout showing the model progression. Warnings will still printed out. + +- **-l --loud output progression given as time step, date and discharge** + The default on-screen output of the lisflood command is the step count and the date/time of each computational step (each step being the run of all the activated modules in sequence for each time step). By setting this "-l" flag, the output will show also the discharge value at each step. + +- **-c --checkfiles input maps and stack maps are checked, output for each input map BUT no model run** + Check the correctness of input maps. The option generates also a table of all the input maps specified in the input xml file. The table will show name of the map, the file path (or single value of the variable), the number of non missing values (nonMV) and missing values (MV) after the comparison with the mask map and LDD map, minimum, maximum and average values for each map. In case any map have missing values, a warning will be issued on screen. The model will NOT run, the lisflood command with -c flags will only check the maps. + +- **-n --nancheck check input maps and routing output for any NaN value generated during model run** + Check for NaN values in input maps and routing output values during the model run. This check do not include missing values (e.g. -9999) in input maps. + +- **-h --noheader .tss file have no header and start immediately with the time series** + By default a header including number of values, value type, settings file path and date of creation is printed when generating the output time series files. Using this flag the header will not be printed in Time Series .tss output files, and the file will start immediately with values. + +- **-d --debug debug outputs** + This flags is used to generate extensive debug txt files containing variable values for each time steps. Each file will be named as "Debug_out_<stepnumber>.txt" + +- **-i --initonly only run initialisation, not the dynamic loop** + By adding this flag the model will only run the initialization and will not enter the dynamic loop (i.e. it stops before computing the first time step) + +- **-s --skipvalreplace skip replacement of invalid values in meteo input maps (ignore valid_min and valid_max)** + By adding this flag the model will skip replacement of invalid values in meteo input maps (ignore valid_min and valid_max). In normal execution, values outside valid_min and valid_max range are replaced with NaN \ No newline at end of file diff --git a/docs/4_annex_tests/index.md b/docs/4_annex_tests/index.md index 55162f34..04343343 100644 --- a/docs/4_annex_tests/index.md +++ b/docs/4_annex_tests/index.md @@ -31,19 +31,23 @@ Some tests use `array_equal` option in order to compare values using [`numpy.arr Tests that are using Comparator classes are: -| Test name | File | Usage and tolerences | -|-------------------------|------------------------|-----------------------------------------------------------------------------------| -| test_dates_steps_day | test_dates_steps.py | NetCDFComparator(array_equal=True) | -| test_dates_steps_6h | test_dates_steps.py | NetCDFComparator(array_equal=True) | -| test_end_state_reported | test_state_end_maps.py | NetCDFComparator(array_equal=True) | -| test_dis_daily | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | -| test_dis_6h | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | -| test_init_daily | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | -| test_init_6h | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | -| test_warmstart_daily | test_warmstart.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(array_equal=True) | -| test_warmstart_6h | test_warmstart.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(array_equal=True) | -| test_subcacthment_daily | test_subcatchments.py | NetCDFComparator(array_equal=True) | -| test_subcacthment_6h | test_subcatchments.py | NetCDFComparator(array_equal=True) | +| Test name | File | Usage and tolerences | +|--------------------------|--------------------------|-----------------------------------------------------------------------------------| +| test_dates_steps_day | test_dates_steps.py | NetCDFComparator(array_equal=True) | +| test_dates_steps_6h | test_dates_steps.py | NetCDFComparator(array_equal=True) | +| test_end_state_reported | test_state_end_maps.py | NetCDFComparator(array_equal=True) | +| test_dis_daily | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | +| test_dis_6h | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | +| test_init_daily | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | +| test_init_6h | test_results.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(atol=0.0001, rtol=0.001) | +| test_warmstart_daily | test_warmstart.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(array_equal=True) | +| test_warmstart_6h | test_warmstart.py | NetCDFComparator(atol=0.0001, rtol=0.001), TSSComparator(array_equal=True) | +| test_subcacthment_daily | test_subcatchments.py | NetCDFComparator(array_equal=True) | +| test_subcacthment_6h | test_subcatchments.py | NetCDFComparator(array_equal=True) | +| test_reported_steps | test_reported_steps.py | NetCDFComparator(array_equal=True) | +| test_waterabstraction_24h| test_water_abstraction.py| NetCDFComparator(array_equal=True) | +| test_waterabstraction_6h | test_water_abstraction.py| NetCDFComparator(array_equal=True) | + ## How to execute tests @@ -316,6 +320,52 @@ def test_rep_dischargetss(self): tss_to_check='disWin.tss') ``` +### Testing reporting steps +Make sure that OSLisflood prints state files following reporting steps formula in ReportSteps xml option. + +#### Implementation +[test_reported_steps.py](https://github.com/ec-jrc/lisflood-code/blob/master/tests/test_reported_steps.py){:target="_blank"} + +In order to test that specific files are written when a step matches the ReportSteps formula, we use a test formula 'starttime+10..endtime'. +Then we geterate specific outputs for the desired steps and compare the two output using NetCDFComparator(array_equal=True) + +**Example: test_reported_steps** + +```python + def test_reported_steps(self): + + strReportStepA = 'starttime+10..endtime' + settings_a = setoptions(self.settings_file, + vars_to_set={'StepStart': '01/07/2016 00:00', 'StepEnd': '01/08/2016 00:00', + 'PathOut': '$(PathRoot)/out/a', + 'ReportSteps': strReportStepA, + }) + mk_path_out(self.out_dir_a) + lisfloodexe(settings_a) + + int_start, _ = datetoint(settings_a.binding['StepStart'], settings_a.binding) + int_end, _ = datetoint(settings_a.binding['StepEnd'], settings_a.binding) + + strReportStepB = '' + for i in range(int_start,int_end,10): + if strReportStepB != '': + strReportStepB = strReportStepB + ',' + strReportStepB = strReportStepB + str(i) + + settings_b = setoptions(self.settings_file, + vars_to_set={'StepStart': '01/07/2016 00:00', 'StepEnd': '01/08/2016 00:00', + 'PathOut': '$(PathRoot)/out/b', + 'ReportSteps': strReportStepB, + }) + mk_path_out(self.out_dir_b) + lisfloodexe(settings_b) + + comparator = NetCDFComparator(settings_a.maskpath, array_equal=True) + comparator.compare_dirs(self.out_dir_b, self.out_dir_a) + assert not comparator.errors +``` + + ### Testing Init run output maps Make sure OSLisflood can run an initial run to generate AVGDIS and LZAVIN maps with proper extension (.nc or .map) @@ -382,6 +432,44 @@ def test_dates_steps_daily(self): comparator.compare_dirs(out_a, out_b) ``` +### Testing Water Abstraction +The test verifies the use of the correct (closest, antecedent timestamp, and closest, antecedent day of the year) water demand information, unsing the following options: +TransientWaterDemandChange = this option allows to read water demand information from a netcdf stack which includes simulation start and end dates. The water demand information for the closest, antecedent timestamp (dd/mm/yyyy) to the modelled time step is used in the computations. +UseWaterDemandAveYear = this option allows to read water demand information from a netcdf stack with a single year temporal coverage. The water demand information for the closest, antecedent day of the year (dd/mm) to the modelled time step is used in the computations. +For more information please refer to [Water use - LISFLOOD (ec-jrc.github.io)](https://ec-jrc.github.io/lisflood-model/2_18_stdLISFLOOD_water-use/) + +#### Implementation +[test_water_abstraction.py](https://github.com/ec-jrc/lisflood-code/blob/master/tests/test_water_abstraction.py){:target="_blank"} +The test uses two datasets: waterdemand19902019, that includes daily data from 1990 to 2019, and waterdemand, that includes only one year used as reference for any year. +Test asserts that output generated using TransientWaterDemandChange with useWaterDemandAveYear flag active using the reference dataset included into the waterdemand folder is the same of the one generated disabling useWaterDemandAveYear flag and using the waterdemand19902019 folder. + + +```python + def run_lisflood_waterabstraction(self, dt_sec): + + settings_a = setoptions(self.settings_file, + opts_to_set=('TransientWaterDemandChange', 'useWaterDemandAveYear'), + vars_to_set={'StepStart': '30/07/2016 00:00', 'StepEnd': '01/08/2016 00:00', + 'DtSec': dt_sec, 'PathOut': '$(PathRoot)/out/a', + 'PathWaterUse': '$(PathRoot)/maps/waterdemand' + }) + mk_path_out(self.out_dir_a) + lisfloodexe(settings_a) + + settings_b = setoptions(self.settings_file, + opts_to_set=('TransientWaterDemandChange'), + opts_to_unset=('useWaterDemandAveYear'), + vars_to_set={'StepStart': '30/07/2016 00:00', 'StepEnd': '01/08/2016 00:00', + 'DtSec': dt_sec, 'PathOut': '$(PathRoot)/out/b', + 'PathWaterUse': '$(PathRoot)/maps/waterdemand19902019'}) + mk_path_out(self.out_dir_b) + lisfloodexe(settings_b) + + comparator = NetCDFComparator(settings_a.maskpath, array_equal=True) + comparator.compare_dirs(self.out_dir_b, self.out_dir_a) + +``` + ## Other LF tests included in repository There are other tests included in [tests/](https://github.com/ec-jrc/lisflood-code/tree/master/tests){:target="_blank"}. diff --git a/docs/_config.yml b/docs/_config.yml index 2afa68ca..806d9cca 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -216,6 +216,8 @@ defaults: url: 3_step6_running-LISFLOOD - title: "Step 6: model output" url: 3_step6_model-output + - title: "Step 7: commmand line flags" + url: 3_step7_commmand-line-flags - section_title: "Static maps" items: diff --git a/docs/media/lisfloodSettings_reference_coldstart_optional_init_v4.0.0.xml b/docs/media/lisfloodSettings_reference_coldstart_optional_init_v4.0.0.xml index ea8dd90d..e8441acd 100644 --- a/docs/media/lisfloodSettings_reference_coldstart_optional_init_v4.0.0.xml +++ b/docs/media/lisfloodSettings_reference_coldstart_optional_init_v4.0.0.xml @@ -99,7 +99,8 @@ - + + @@ -133,23 +134,22 @@ The option "MapsCaching" may take the following values: -************************************************************** -PARALLELIZED KINEMATIC ROUTING -************************************************************** + +The option "OutputMapsChunks" may take the following values: + - "[positive integer number]" : Dump outputs to disk every X steps (default 1) + -!-- Parallelisation of the kinematic wave routing (via openMP multi-threading). -The option "numCPUs_parallelKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) - (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> - + +The option "OutputMapsDataType" sets the output data type and may take the following values: + - "float64" + - "float32" + + ************************************************************** -PARALLELISATION WITH NUMBA (USED IN SOILLOOP) +PARALLELISATION WITH NUMBA (USED IN ROUTING AND SOILLOOP) ************************************************************** @@ -1692,7 +1692,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT ************************************************************** --> - @@ -1702,6 +1701,11 @@ This is necessary when using projected coordinates (x,y) + + + + + Clone map diff --git a/docs/media/lisfloodSettings_reference_prerun_optional_init_v4.0.0.xml b/docs/media/lisfloodSettings_reference_prerun_optional_init_v4.0.0.xml index 96dfdaea..ae92afed 100644 --- a/docs/media/lisfloodSettings_reference_prerun_optional_init_v4.0.0.xml +++ b/docs/media/lisfloodSettings_reference_prerun_optional_init_v4.0.0.xml @@ -80,7 +80,7 @@ # report end maps - + # report maps @@ -99,7 +99,8 @@ - + + @@ -133,23 +134,22 @@ The option "MapsCaching" may take the following values: -************************************************************** -PARALLELIZED KINEMATIC ROUTING -************************************************************** + +The option "OutputMapsChunks" may take the following values: + - "[positive integer number]" : Dump outputs to disk every X steps (default 1) + -!-- Parallelisation of the kinematic wave routing (via openMP multi-threading). -The option "numCPUs_parallelKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) - (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> - + +The option "OutputMapsDataType" sets the output data type and may take the following values: + - "float64" + - "float32" + + ************************************************************** -PARALLELISATION WITH NUMBA (USED IN SOILLOOP) +PARALLELISATION WITH NUMBA (USED IN ROUTING AND SOILLOOP) ************************************************************** @@ -1692,7 +1692,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT ************************************************************** --> - @@ -1702,6 +1701,11 @@ This is necessary when using projected coordinates (x,y) + + + + + Clone map diff --git a/environment.yml b/environment.yml index 20a42945..571652d8 100644 --- a/environment.yml +++ b/environment.yml @@ -149,7 +149,6 @@ dependencies: - cftime==1.6.0 - cloudpickle==2.1.0 - coverage==6.4.1 - - cython==0.29.30 - dask==2022.2.0 - fsspec==2022.5.0 - future==0.18.2 diff --git a/requirements.txt b/requirements.txt index f1dc6e0e..297cc0c0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,6 @@ attrs>=19.3.0 beautifulsoup4>=4.8.1 cftime>=1.0.4.2 coverage>=6.0 -Cython>=0.29.14 dask>=2021.10.0 future>=0.18.0 GDAL>=3.2.0 @@ -15,7 +14,7 @@ more-itertools>=7.2.0 netCDF4>=1.5.3,<=1.6.4 nine>=1.0.0 numba>=0.54.0 -numexpr==2.8.1 +numexpr>=2.7.0,<=2.8.4 numpy>=1.21.0 packaging>=19.2 pandas>=1,<2 diff --git a/setup.py b/setup.py index 985e4d99..f69f1e17 100755 --- a/setup.py +++ b/setup.py @@ -59,27 +59,6 @@ version = f.read().strip() print(">>>>>>>>>>>>>>>> Building LISFLOOD version {} <<<<<<<<<<<<<<<<<<".format(version)) -try: - # noqa - from Cython.Build import cythonize - HAS_CYTHON = True -except ImportError: - HAS_CYTHON = False - -src_ext = 'src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.{}' -extension_ext = 'pyx' if HAS_CYTHON else 'c' - -ext_modules = [ - Extension( - 'lisflood.hydrological_modules.kinematic_wave_parallel_tools', - [src_ext.format(extension_ext)], - extra_compile_args=['-O3', "-march=native", '-fopenmp'], - extra_link_args=['-fopenmp'] - ) -] - -if HAS_CYTHON: - ext_modules = cythonize(ext_modules) readme_file = os.path.join(current_dir, 'README.md') with open(readme_file, 'r') as f: @@ -177,11 +156,9 @@ def _get_gdal_version(): 'nine', 'setuptools>=41.0', 'numpy>=1.15', - 'cython', ], install_requires=requirements, scripts=['bin/lisflood'], - ext_modules=ext_modules, zip_safe=True, classifiers=[ # complete classifier list: http://pypi.python.org/pypi?%3Aaction=list_classifiers diff --git a/src/lisflood/Lisflood_dynamic.py b/src/lisflood/Lisflood_dynamic.py index c76673fd..70f4c188 100644 --- a/src/lisflood/Lisflood_dynamic.py +++ b/src/lisflood/Lisflood_dynamic.py @@ -56,11 +56,12 @@ def dynamic(self): self.TimeSinceStart = self.currentTimeStep() - self.firstTimeStep() + 1 if flags['loud']: - print("%-6i %10s" % (self.currentTimeStep(), self.CalendarDate.strftime("%d/%m/%Y %H:%M"))) + sys.stdout.write("%-6i %20s" % (self.currentTimeStep(), self.CalendarDate.strftime("%d/%m/%Y %H:%M"))) else: if not flags['checkfiles']: if flags['quiet'] and not flags['veryquiet']: sys.stdout.write(".") + sys.stdout.flush() if not flags['quiet'] and not flags['veryquiet']: # Print step number and date to console sys.stdout.write("\r%d" % i), sys.stdout.write("%s" % " - "+self.CalendarDate.strftime("%d/%m/%Y %H:%M")) diff --git a/src/lisflood/__init__.py b/src/lisflood/__init__.py index dd927604..524386a0 100644 --- a/src/lisflood/__init__.py +++ b/src/lisflood/__init__.py @@ -10,8 +10,8 @@ __version__ = version __authors__ = "Ad de Roo, Emiliano Gelati, Peter Burek, Johan van der Knijff, Niko Wanders" -__date__ = "15/11/2023" -__copyright__ = "Copyright 2019-2023, European Commission - Joint Research Centre" +__date__ = "05/01/2024" +__copyright__ = "Copyright 2019-2024, European Commission - Joint Research Centre" __maintainer__ = "Stefania Grimaldi, Cinzia Mazzetti, Carlo Russo, Damien Decremer, Corentin Carton De Wiart, Valerio Lorini, Ad de Roo" __status__ = "Operation" __institution__ = "European Commission DG Joint Research Centre (JRC) - E1, D2 Units" diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py index a84b9963..e787c5ac 100644 --- a/src/lisflood/global_modules/default_options.py +++ b/src/lisflood/global_modules/default_options.py @@ -502,7 +502,7 @@ 'LakePrevInflowEnd': ReportedMap(name='LakePrevInflowEnd', output_var='LakeInflowOld', unit='m3/s', end=['repEndMaps'], steps=[], all=[], - restrictoption=['simulateLakes'], + restrictoption=['nonInit', 'simulateLakes'], monthly=False, yearly=False), 'LakePrevInflowState': ReportedMap(name='LakePrevInflowState', output_var='LakeInflowOld', unit='m3/s', end=[], steps=['repStateMaps'], @@ -511,7 +511,7 @@ monthly=False, yearly=False), 'LakePrevOutflowEnd': ReportedMap(name='LakePrevOutflowEnd', output_var='LakeOutflow', unit='m3/s', end=['repEndMaps'], steps=[], all=[], - restrictoption=['simulateLakes'], + restrictoption=['nonInit', 'simulateLakes'], monthly=False, yearly=False), 'LakePrevOutflowState': ReportedMap(name='LakePrevOutflowState', output_var='LakeOutflow', unit='m3/s', end=[], steps=['repStateMaps'], diff --git a/src/lisflood/global_modules/output.py b/src/lisflood/global_modules/output.py index 11027027..aeb8acf3 100644 --- a/src/lisflood/global_modules/output.py +++ b/src/lisflood/global_modules/output.py @@ -17,6 +17,7 @@ import os import numpy as np from pcraster import ifthen, catchmenttotal, mapmaximum +import sys from .zusatz import TimeoutputTimeseries from .add1 import decompress, valuecell, loadmap, compressArray @@ -556,9 +557,10 @@ def dynamic(self): if flags['loud']: # print the discharge of the first output map loc try: - print(" %10.2f" % self.var.Tss["DisTS"].firstout(decompress(self.var.ChanQAvg))) + sys.stdout.write(" %10.2f" % self.var.Tss["DisTS"].firstout(decompress(self.var.ChanQAvg))) except: pass + sys.stdout.write("\n") for tss in report_time_serie_act: # report time series diff --git a/src/lisflood/hydrological_modules/compile_kinematic_wave_parallel_tools.py b/src/lisflood/hydrological_modules/compile_kinematic_wave_parallel_tools.py deleted file mode 100755 index f94ee4f5..00000000 --- a/src/lisflood/hydrological_modules/compile_kinematic_wave_parallel_tools.py +++ /dev/null @@ -1,34 +0,0 @@ -""" - -Copyright 2019 European Union - -Licensed under the EUPL, Version 1.2 or as soon they will be approved by the European Commission subsequent versions of the EUPL (the "Licence"); - -You may not use this work except in compliance with the Licence. -You may obtain a copy of the Licence at: - -https://joinup.ec.europa.eu/sites/default/files/inline-files/EUPL%20v1_2%20EN(1).txt - -Unless required by applicable law or agreed to in writing, software distributed under the Licence is distributed on an "AS IS" basis, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the Licence for the specific language governing permissions and limitations under the Licence. - -""" - -from distutils.core import setup -from distutils.extension import Extension -from Cython.Build import cythonize - -ext_modules = [ - Extension( - "kinematic_wave_parallel_tools", - ["kinematic_wave_parallel_tools.pyx"], - extra_compile_args=["-O3", "-march=native", "-fopenmp"], - extra_link_args=["-fopenmp"] - ) -] - -setup( - name='kinematic_wave-parallel', - ext_modules=cythonize(ext_modules, annotate=True), -) diff --git a/src/lisflood/hydrological_modules/compile_kinematic_wave_tools.md b/src/lisflood/hydrological_modules/compile_kinematic_wave_tools.md deleted file mode 100755 index 9ba572b4..00000000 --- a/src/lisflood/hydrological_modules/compile_kinematic_wave_tools.md +++ /dev/null @@ -1,20 +0,0 @@ -# Lisflood code - -Source code of lisflood model. - -To compile this Cython module to enable OpenMP multithreading (parallel kinematic wave): - -1) Delete the files *.so (if any) in directory hydrological-modules - -2) Inside the hydrological_modules folder, execute "python compile_kinematic_wave_parallel_tools.py build_ext --inplace" - -Important: the module has to be compiled on the machine where the model is run - the resulting binary is not portable. - -Then in the settings file the option "numberParallelThreadsKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> -```xml - -``` diff --git a/src/lisflood/hydrological_modules/kinematic_wave_parallel.py b/src/lisflood/hydrological_modules/kinematic_wave_parallel.py index 2fe25659..48f51cfd 100755 --- a/src/lisflood/hydrological_modules/kinematic_wave_parallel.py +++ b/src/lisflood/hydrological_modules/kinematic_wave_parallel.py @@ -15,8 +15,8 @@ Parallel implementation of the kinematic wave routing. This module provides tools to pre-process the flow direction matrix, and the kinematicWave class that -solves the kinematic wave equation. Most functions and methods call functions in the Cython module -kinematic_wave_parallel_tools.pyx, for performance reasons. +solves the kinematic wave equation. Most functions and methods call functions in the Numba module +kinematic_wave_parallel_tools.py, for performance reasons. Parallelisation is achieved by grouping the pixels in ordered sets. Within each set, provided that pixelsm in previous sets have already been routed, pixels can be routed independently of each other and thus in parallel. For further details, see the method kinematicWave._setRoutingOrders. @@ -40,39 +40,8 @@ from ..global_modules.errors import LisfloodWarning -# IMPORT THE PARALLELISE KINEMATIC WAVE RUTING MODULE: IF IT WAS NOT COMPILED ON THE CURRENT MACHINE, -# ROUTING IS EXECUTED SERIALLY if the binary .so file does not exist or is not newer than the source .pyx, -# then the Cython module is imported directly from the source. -# For safety against binary files compiled on other machines and copied here, -# the binary must be at least 10 seconds younger than the source (see line 13). -# Importing directly from the source prevents using OpenMP multithreading. -# In such case, the routing is executed serially. - -WINDOWS_OS = system() == "Windows" -ROOT = os.path.join(os.path.dirname(os.path.realpath(__file__)), "kinematic_wave_parallel_tools") -SRC = '{}.pyx'.format(ROOT) -BINS = '{}*.so'.format(ROOT) if not WINDOWS_OS else '{}*.pyd'.format(ROOT) - -bins = glob.glob(BINS) - -older = False -for binary in bins: - if os.stat(binary).st_mtime < os.stat(SRC).st_mtime: - older = True - break - -if not bins or older: - # Activate the direct import from source of Cython modules. - import pyximport - setup_args = {"script_args": ["--compiler=mingw32"]} if WINDOWS_OS else None - # If this is executed, the binary .so file will be ignored and the routing is executed serially. - pyximport.install(setup_args=setup_args) - warnings.warn(LisfloodWarning("""The Cython module {} has not been compiled on the current machine (to compile, see instructions in the module's docstring). -The kinematic wave routing is executed serially (not in parallel).""".format(SRC))) - from . import kinematic_wave_parallel_tools as kwpt - # ------------------------------------------------------------------------------------------------- # CONSTANTS # ------------------------------------------------------------------------------------------------- @@ -118,7 +87,7 @@ def streamLookups(flow_dir, land_mask): land_points[land_mask] = np.arange(num_pixs, dtype=int) downstream_lookup, upstream_lookup = kwpt.upDownLookups(flow_dir, np.ascontiguousarray(land_mask).astype(np.uint8), land_points, num_pixs, IX_ADDS) max_num_ups_pixs = max(1, np.any(upstream_lookup != -1, 0).sum()) # maximum number of upstreams pixels - return downstream_lookup, np.ascontiguousarray(upstream_lookup[:,:max_num_ups_pixs]).astype(int) # astype for cython import in windows (see below) + return downstream_lookup, np.ascontiguousarray(upstream_lookup[:,:max_num_ups_pixs]).astype(int) def topoDistFromSea(downstream_lookup, upstream_lookup): """""" @@ -145,8 +114,11 @@ def topoDistFromSea(downstream_lookup, upstream_lookup): class kinematicWave: """""" - def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space_delta, time_delta, num_threads, alpha_floodplains=None): + def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space_delta, time_delta, alpha_floodplains=None, flagnancheck=False): """""" + # variable to avoid printing repeated warning messages + self.kinematic_wave_warning_printed=False + self.flagnancheck=flagnancheck # Parameters for the solution of the discretised Kinematic wave continuity equation self.space_delta = space_delta self.beta = beta @@ -154,8 +126,6 @@ def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space self.b_minus_1 = beta - 1 self.a_dx_div_dt_channel = alpha_channel * space_delta / time_delta self.b_a_dx_div_dt_channel = beta * self.a_dx_div_dt_channel - # Set number of parallel threads (openMP) - self.num_threads = int(num_threads) if 0 < num_threads <= cpu_count() else cpu_count() # If split-routing (floodplains) if alpha_floodplains is not None: self.a_dx_div_dt_floodplains = alpha_floodplains * space_delta / time_delta @@ -163,7 +133,7 @@ def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space # Process flow direction matrix: downstream and upstream lookups, and routing orders flow_dir = decodeFlowMatrix(rebuildFlowMatrix(compressed_encoded_ldd, land_mask)) self.downstream_lookup, self.upstream_lookup = streamLookups(flow_dir, land_mask) - self.num_upstream_pixels = (self.upstream_lookup != -1).sum(1).astype(int) # astype for cython import in windows (to avoid 'long long' buffer dtype mismatch) + self.num_upstream_pixels = (self.upstream_lookup != -1).sum(1).astype(int) # Routing order: decompose domain into batches; within each batch, pixels can be routed in parallel self._setRoutingOrders() @@ -184,8 +154,8 @@ def _setRoutingOrders(self): self.pixels_ordered = self.pixels_ordered.sort(["order", "pixels"]).set_index("order").squeeze() order_counts = self.pixels_ordered.groupby(self.pixels_ordered.index).count() stop = order_counts.cumsum() - self.order_start_stop = np.column_stack((np.append(0, stop[:-1]), stop)).astype(int) # astype for cython import in windows (see above) - self.pixels_ordered = self.pixels_ordered.values.astype(int) # astype for cython import in windows (see above) + self.order_start_stop = np.column_stack((np.append(0, stop[:-1]), stop)).astype(int) + self.pixels_ordered = self.pixels_ordered.values.astype(int) def kinematicWaveRouting(self, discharge, specific_lateral_inflow, section="main_channel"): """Kinematic wave routing: wrapper around kinematic_wave_parallel_tools.kinematicWave""" @@ -206,4 +176,9 @@ def kinematicWaveRouting(self, discharge, specific_lateral_inflow, section="main # Solve the Kinematic wave equation kwpt.kinematicRouting(discharge, lateral_inflow, constant, self.upstream_lookup,\ self.num_upstream_pixels, self.pixels_ordered, self.order_start_stop,\ - self.beta, self.inv_beta, self.b_minus_1, a_dx_div_dt, b_a_dx_div_dt, self.num_threads) + self.beta, self.inv_beta, self.b_minus_1, a_dx_div_dt, b_a_dx_div_dt) + if self.flagnancheck: + if self.kinematic_wave_warning_printed==False: + if np.all(np.isfinite(discharge))==False: + warnings.warn(LisfloodWarning("Warning: NaN or Inf values after kinematicRouting module. Suggestion: please check the input maps (e.g. channel geometry and ldd)")) + self.kinematic_wave_warning_printed=True \ No newline at end of file diff --git a/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.pyx b/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py similarity index 60% rename from src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.pyx rename to src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py index d3ebeafa..8646cefa 100755 --- a/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.pyx +++ b/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py @@ -1,4 +1,3 @@ -#cython: language_level=2, boundscheck=False """ Copyright 2019 European Union @@ -14,58 +13,46 @@ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the Licence for the specific language governing permissions and limitations under the Licence. -Collection of functions called by the Python module kinematic_wave_parallel.py. These functions are -implemented in Cython to achieve compiled-code performance when executing loops over numpy arrays. -To compile this Cython module to enable OpenMP multithreading, execute: -$ python compile_kinematic_wave_parallel_tools.py build_ext --inplace -Important: the module has to be compiled on the machine where the model is run - the resulting binary is not portable. +Collection of functions called by the Python module kinematic_wave_parallel.py. """ -__author__ = "Emiliano Gelati" -__contact__ = "emiliano.gelati@ec.europa.eu" -__version__ = "1.0" -__date__ = "2016/06/01" - +from math import fabs import numpy as np -from cython import boundscheck, wraparound, cdivision -from cython.parallel import prange -from libc.math cimport fabs, fmax, fmin +from numba import njit, prange +from builtins import max + + -cdef: - double NEWTON_TOL = 1e-12 # Max allowed error in iterative solution - long MAX_ITERS = 3000 # Max number of iterations - double MIN_DISCHARGE = 1e-30 # Minimum allowed discharge value to avoid numerical instability +NEWTON_TOL = 1e-12 # Max allowed error in iterative solution +MAX_ITERS = 3000 # Max number of iterations +MIN_DISCHARGE = 1e-30 # Minimum allowed discharge value to avoid numerical instability # ------------------------------------------------------------------------------------------------- # ROUTING FUNCTIONS # ------------------------------------------------------------------------------------------------- -@boundscheck(False) -@wraparound(False) -def kinematicRouting(double[::1] discharge, double[::1] lateral_inflow, double[::1] constant, long[:,::1] upstream_lookup,\ - long[::1] num_upstream_pixels, long[::1] ordered_pixels, long[:,::1] start_stop, double beta, double inv_beta,\ - double b_minus_1, double[::1] a_dx_div_dt, double[::1] b_a_dx_div_dt, long num_threads): +@njit(parallel=True, fastmath=False, cache=True) +def kinematicRouting(discharge, lateral_inflow, constant, upstream_lookup,\ + num_upstream_pixels, ordered_pixels, start_stop, beta, inv_beta,\ + b_minus_1, a_dx_div_dt, b_a_dx_div_dt): """""" - cdef: - long order, index, first, last, num_orders = start_stop.shape[0] + num_orders = start_stop.shape[0] # Iterate through routing orders (sets of pixels for which the kinemativc wave can be solved independently and thus in parallel) for order in range(num_orders): first = start_stop[order,0] last = start_stop[order,1] - for index in prange(first, last, nogil=True, schedule="dynamic", num_threads=num_threads): + for index in prange(first, last): solve1Pixel(ordered_pixels[index], discharge, lateral_inflow, constant, upstream_lookup,\ num_upstream_pixels, a_dx_div_dt, b_a_dx_div_dt, beta, inv_beta, b_minus_1) -@boundscheck(False) -@wraparound(False) -@cdivision(True) -cdef void solve1Pixel(long pix, double[::1] discharge, double[::1] lateral_inflow, double[::1] constant,\ - long[:,::1] upstream_lookup, long[::1] num_upstream_pixels, double[::1] a_dx_div_dt,\ - double[::1] b_a_dx_div_dt, double beta, double inv_beta, double b_minus_1) nogil: +@njit(nogil=True, fastmath=False, cache=True) +def solve1Pixel(pix, discharge, lateral_inflow, constant,\ + upstream_lookup, num_upstream_pixels, a_dx_div_dt,\ + b_a_dx_div_dt, beta, inv_beta, b_minus_1): """""" - cdef: - long ups_ix, count = 0 - double error, const_plus_ups_infl, a_cpui_pow_b_m_1, secant_bound, other_bound, previous_estimate = -1., upstream_inflow = 0. + count = 0 + previous_estimate = -1.0 + upstream_inflow = 0.0 # Inflow from upstream pixels for ups_ix in range(num_upstream_pixels[pix]): upstream_inflow += discharge[upstream_lookup[pix,ups_ix]] @@ -87,16 +74,20 @@ def kinematicRouting(double[::1] discharge, double[::1] lateral_inflow, double[: while fabs(error) > NEWTON_TOL and discharge[pix] != previous_estimate and count < MAX_ITERS: # is previous_estimate useful? previous_estimate = discharge[pix] discharge[pix] -= error / (1 + b_a_dx_div_dt[pix] * discharge[pix]**b_minus_1) - discharge[pix] = fmax(discharge[pix], NEWTON_TOL) + discharge[pix] = max(discharge[pix], NEWTON_TOL) error = closureError(discharge[pix], const_plus_ups_infl, a_dx_div_dt[pix], beta) count += 1 # If iterations converge to NEWTON_TOL, set value to 0 if discharge[pix] == NEWTON_TOL: discharge[pix] = 0 - -@boundscheck(False) -@wraparound(False) -cdef double closureError(double discharge, double upper_bound, double a_dx_div_dt, double beta) nogil: + # to simulate inf or nan: discharge[pix] = 1.0/0.0 + # with gil: + # got_valid_value = np.isfinite(discharge[pix]) + # if got_valid_value==False: + # print(str(discharge[pix]) + ' found at pix:' + str(pix)) + +@njit(nogil=True, fastmath=False, cache=True) +def closureError(discharge, upper_bound, a_dx_div_dt, beta): """""" return discharge + a_dx_div_dt * discharge**beta - upper_bound @@ -105,11 +96,10 @@ def kinematicRouting(double[::1] discharge, double[::1] lateral_inflow, double[: # ------------------------------------------------------------------------------------------------- # FLOW DIRECTION MATRIX PRE-PROCESSING FUNCTIONS # ------------------------------------------------------------------------------------------------- -def updateUpstreamRouted(long[::1] pixels_routed, long[::1] downstream_lookup,\ - long[:,::1] upstream_lookup, unsigned char[:,::1] upstream_routed): +@njit(nogil=True, fastmath=False, cache=True) +def updateUpstreamRouted(pixels_routed, downstream_lookup,\ + upstream_lookup, upstream_routed): """Called by kinematic_wave_parallel.orderKinematicRouting to implement greedy pixel ordering starting from headwaters.""" - cdef: - long downs_pix, pix, ups_index for pix in pixels_routed: downs_pix = downstream_lookup[pix] if downs_pix != -1: @@ -118,14 +108,14 @@ def updateUpstreamRouted(long[::1] pixels_routed, long[::1] downstream_lookup,\ ups_index += 1 upstream_routed[downs_pix,ups_index] = True -def upDownLookups(long[:,::1] flow_d8, unsigned char[:,::1] land_mask, long[:,::1] land_points, long num_pixs, long[:,::1] ix_adds): +@njit(nogil=True, fastmath=False, cache=True) +def upDownLookups(flow_d8, land_mask, land_points, num_pixs, ix_adds): '''Called by catchment.streamLookups''' - cdef: - long[::1] downstream_lookup = -np.ones(num_pixs, int) - long[:,::1] upstream_lookup = -np.ones((num_pixs, 8), int) - unsigned int[::1] ups_count = np.zeros(num_pixs, np.uintc) - int src_row, src_col, dst_row, dst_col - unsigned int ups_p, downs_p, num_row = flow_d8.shape[0], num_col = flow_d8.shape[1] + downstream_lookup = -np.ones(num_pixs) + upstream_lookup = -np.ones((num_pixs, 8)) + ups_count = np.zeros(num_pixs, np.uintc) + num_row = flow_d8.shape[0] + num_col = flow_d8.shape[1] for src_row in range(num_row): for src_col in range(num_col): if flow_d8[src_row,src_col] < 8: @@ -144,12 +134,13 @@ def upDownLookups(long[:,::1] flow_d8, unsigned char[:,::1] land_mask, long[:,:: # ------------------------------------------------------------------------------------------------- # MISCELLANOUS TOOLS # ------------------------------------------------------------------------------------------------- -def immediateUpstreamInflow(double[::1] discharge, long[:,::1] upstream_lookup, long[::1] num_upstream_pixels): +@njit(nogil=True, fastmath=False, cache=True) +def immediateUpstreamInflow(discharge, upstream_lookup, num_upstream_pixels): """""" - cdef: - long pix, ups_ix, num_pixels = discharge.size - double [::1] inflow = np.zeros(num_pixels) + num_pixels = discharge.size + inflow = np.zeros(num_pixels) for pix in range(num_pixels): for ups_ix in range(num_upstream_pixels[pix]): inflow[pix] += discharge[upstream_lookup[pix,ups_ix]] return np.asarray(inflow) + diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py index 4683437e..aa64c110 100644 --- a/src/lisflood/hydrological_modules/routing.py +++ b/src/lisflood/hydrological_modules/routing.py @@ -343,7 +343,8 @@ def initialSecond(self): """ settings = LisSettings.instance() option = settings.options - binding = settings.binding + flags = settings.flags + self.var.ChannelAlpha2 = None # default value, if split-routing is not active and only water is routed only in the main channel # ************************************************************ @@ -399,7 +400,7 @@ def initialSecond(self): maskinfo = MaskInfo.instance() self.river_router = kinematicWave(compressArray(self.var.LddKinematic), ~maskinfo.info.mask, self.var.ChannelAlpha, self.var.Beta, self.var.ChanLength, self.var.DtRouting, - int(binding["numCPUs_parallelKinematicWave"]), alpha_floodplains=self.var.ChannelAlpha2) + alpha_floodplains=self.var.ChannelAlpha2, flagnancheck=flags['nancheck']) if option['InitLisflood'] and option['repMBTs']: self.var.StorageStepINIT= self.var.ChanM3Kin @@ -436,7 +437,6 @@ def dynamic(self, NoRoutingExecuted): """ settings = LisSettings.instance() option = settings.options - binding = settings.binding if not(option['InitLisflood']): # only with no InitLisflood self.lakes_module.dynamic_inloop(NoRoutingExecuted) diff --git a/src/lisflood/hydrological_modules/surface_routing.py b/src/lisflood/hydrological_modules/surface_routing.py index 200c8425..00d54faa 100644 --- a/src/lisflood/hydrological_modules/surface_routing.py +++ b/src/lisflood/hydrological_modules/surface_routing.py @@ -103,15 +103,14 @@ def initialSecond(self): maskinfo = MaskInfo.instance() land_mask = ~maskinfo.info.mask settings = LisSettings.instance() - binding = settings.binding - num_threads = int(binding["numCPUs_parallelKinematicWave"]) + flags = settings.flags self.direct_surface_router = kinematicWave(compressArray(self.var.LddToChan), land_mask, self.var.OFAlpha.values[self.var.dim_runoff[1].index('Direct')], self.var.Beta,\ - self.var.PixelLength, dt_surf_routing, num_threads) + self.var.PixelLength, dt_surf_routing, flagnancheck=flags['nancheck']) self.other_surface_router = kinematicWave(compressArray(self.var.LddToChan), land_mask, self.var.OFAlpha.values[self.var.dim_runoff[1].index('Other')], self.var.Beta,\ - self.var.PixelLength, dt_surf_routing, num_threads) + self.var.PixelLength, dt_surf_routing, flagnancheck=flags['nancheck']) self.forest_surface_router = kinematicWave(compressArray(self.var.LddToChan), land_mask, self.var.OFAlpha.values[self.var.dim_runoff[1].index('Forest')], self.var.Beta,\ - self.var.PixelLength, dt_surf_routing, num_threads) + self.var.PixelLength, dt_surf_routing, flagnancheck=flags['nancheck']) def dynamic(self): """ dynamic part of the surface routing module diff --git a/src/lisflood/main.py b/src/lisflood/main.py index bfff7425..034e78bf 100755 --- a/src/lisflood/main.py +++ b/src/lisflood/main.py @@ -85,31 +85,6 @@ def lisfloodexe(lissettings=None): for key in report_steps: report_steps[key] = [x for x in report_steps[key] if model_steps[0] <= x <= model_steps[1]] - if not flags['veryquiet']: - print("\nStart Step - End Step: ", model_steps[0], " - ", model_steps[1]) - print("Start Date - End Date: ", - inttodate(model_steps[0] - 1, calendar(binding['CalendarDayStart'], binding['calendar_type'])), - " - ", - inttodate(model_steps[1] - 1, calendar(binding['CalendarDayStart'], binding['calendar_type']))) - - if flags['loud']: - # print state file date - print("State file Date: ") - try: - print(inttodate(calendar(binding["timestepInit"], binding['calendar_type']), - calendar(binding['CalendarDayStart'], binding['calendar_type']))) - except: - print(calendar(binding["timestepInit"], binding['calendar_type'])) - - # CM: print start step and end step for reporting model state maps - print("Start Rep Step - End Rep Step: ", report_steps['rep'][0], " - ", report_steps['rep'][-1]) - print("Start Rep Date - End Rep Date: ", - inttodate(calendar(report_steps['rep'][0] - 1, binding['calendar_type']), calendar(binding['CalendarDayStart'], binding['calendar_type'])), - " - ", - inttodate(calendar(report_steps['rep'][-1] - 1), calendar(binding['CalendarDayStart'], binding['calendar_type']))) - # messages at model start - print("%-6s %10s %11s\n" % ("Step", "Date", "Discharge")) - # Lisflood is an instance of the class LisfloodModel # LisfloodModel includes 2 methods : initial and dynamic (formulas applied at every timestep) Lisflood = LisfloodModel() @@ -119,6 +94,7 @@ def lisfloodexe(lissettings=None): stLisflood.rquiet = True stLisflood.rtrace = False + if lissettings.mc_set: """ ---------------------------------------------- @@ -148,6 +124,32 @@ def lisfloodexe(lissettings=None): # print info about execution if not flags['veryquiet']: print(LisfloodRunInfo(model_to_run)) + + if not flags['veryquiet']: + print("\nStart Step - End Step: ", model_steps[0], " - ", model_steps[1]) + print("Start Date - End Date: ", + inttodate(model_steps[0] - 1, calendar(binding['CalendarDayStart'], binding['calendar_type'])), + " - ", + inttodate(model_steps[1] - 1, calendar(binding['CalendarDayStart'], binding['calendar_type']))) + + if flags['loud']: + # print state file date + print("State file Date: ") + try: + print(inttodate(calendar(binding["timestepInit"], binding['calendar_type']), + calendar(binding['CalendarDayStart'], binding['calendar_type']))) + except: + print(calendar(binding["timestepInit"], binding['calendar_type'])) + + # CM: print start step and end step for reporting model state maps + print("Start Rep Step - End Rep Step: ", report_steps['rep'][0], " - ", report_steps['rep'][-1]) + print("Start Rep Date - End Rep Date: ", + inttodate(calendar(report_steps['rep'][0] - 1, binding['calendar_type']), calendar(binding['CalendarDayStart'], binding['calendar_type'])), + " - ", + inttodate(calendar(report_steps['rep'][-1] - 1), calendar(binding['CalendarDayStart'], binding['calendar_type']))) + # messages at model start + print("%-6s %15s %18s" % ("Step", "Date", "Discharge")) + # actual run of the model if flags['initonly']: print('initonly flag activated... Stopping now before entering time loop.') @@ -177,7 +179,7 @@ def usage(): -v --veryquiet no output progression is given -l --loud output progression given as time step, date and discharge -c --checkfiles input maps and stack maps are checked, output for each input map BUT no model run - -n --nancheck check NaN values in output maps + -n --nancheck check input maps and routing output for any NaN value generated during model run -h --noheader .tss file have no header and start immediately with the time series -d --debug debug outputs -i --initonly only run initialisation, not the dynamic loop diff --git a/src/lisfloodSettings_reference.xml b/src/lisfloodSettings_reference.xml index 94ea87db..15aceb45 100644 --- a/src/lisfloodSettings_reference.xml +++ b/src/lisfloodSettings_reference.xml @@ -108,7 +108,8 @@ You can use builtin path variables in this template and reference to other paths - + + @@ -156,22 +157,7 @@ The option "OutputMapsDataType" sets the output data type and may take the follo ************************************************************** -PARALLELIZED KINEMATIC ROUTING -************************************************************** - - -!-- Parallelisation of the kinematic wave routing (via openMP multi-threading). -The option "numCPUs_parallelKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) - (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> - - - -************************************************************** -PARALLELISATION WITH NUMBA (USED IN SOILLOOP) +PARALLELISATION WITH NUMBA (USED IN ROUTING AND SOILLOOP) ************************************************************** @@ -1714,7 +1700,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT ************************************************************** --> - diff --git a/tests/data/LF_ETRS89_UseCase/settings/base.xml b/tests/data/LF_ETRS89_UseCase/settings/base.xml index 4c4a14ab..a8faedb7 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/base.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/base.xml @@ -105,15 +105,6 @@ - - - - diff --git a/tests/data/LF_ETRS89_UseCase/settings/cold.xml b/tests/data/LF_ETRS89_UseCase/settings/cold.xml index 06a42eb7..9f789c8c 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/cold.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/cold.xml @@ -142,22 +142,7 @@ ************************************************************** - PARALLELIZED KINEMATIC ROUTING - ************************************************************** - - - !-- Parallelisation of the kinematic wave routing (via openMP multi-threading). - The option "numberParallelThreadsKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) - (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> - - - - ************************************************************** - PARALLELISATION WITH NUMBA (USED IN SOILLOOP) + PARALLELISATION WITH NUMBA (USED IN ROUTING AND SOILLOOP) ************************************************************** @@ -1588,7 +1573,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT ************************************************************** --> - diff --git a/tests/data/LF_ETRS89_UseCase/settings/full.xml b/tests/data/LF_ETRS89_UseCase/settings/full.xml index ad6933e8..66ce403e 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/full.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/full.xml @@ -88,11 +88,6 @@ - "[positive integer number]" : Dump outputs to disk every X steps (default 1) - - ************************************************************** - PARALLELIZED KINEMATIC ROUTING - ************************************************************** - The option "OutputMapsDataType" sets the output data type and may take the following values: @@ -101,18 +96,9 @@ - - - ************************************************************** - PARALLELISATION WITH NUMBA (USED IN SOILLOOP) + PARALLELISATION WITH NUMBA ************************************************************** @@ -2317,7 +2303,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT - diff --git a/tests/data/LF_ETRS89_UseCase/settings/inflow.xml b/tests/data/LF_ETRS89_UseCase/settings/inflow.xml index d94a8aab..f67edfc9 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/inflow.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/inflow.xml @@ -103,15 +103,6 @@ - - - - diff --git a/tests/data/LF_ETRS89_UseCase/settings/prerun.xml b/tests/data/LF_ETRS89_UseCase/settings/prerun.xml index 591d7d16..3dcb0af0 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/prerun.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/prerun.xml @@ -101,15 +101,6 @@ - - - - diff --git a/tests/data/LF_ETRS89_UseCase/settings/warm.xml b/tests/data/LF_ETRS89_UseCase/settings/warm.xml index d3b2ae0c..2fbc9b86 100644 --- a/tests/data/LF_ETRS89_UseCase/settings/warm.xml +++ b/tests/data/LF_ETRS89_UseCase/settings/warm.xml @@ -128,21 +128,6 @@ - - ************************************************************** - PARALLELIZED KINEMATIC ROUTING - ************************************************************** - - - !-- Parallelisation of the kinematic wave routing (via openMP multi-threading). - The option "numberParallelThreadsKinematicWave" may take the following values: - - "0" : auto-detection of the machine/node's number of CPUs (all CPUs are used minus 1) - (do not set it if other simulations are running on the same machine/node) - - "1" : serial execution (not parallel) - - "2", "3", ... : manual setting of the number of parallel threads. - (if exceeding the number of CPUs, the option is set to "0") --> - - ************************************************************** PARALLELISATION WITH NUMBA (USED IN SOILLOOP) @@ -1565,7 +1550,6 @@ LFBINDING: MORE LOW-LEVEL CONTROL OVER MODEL IN- AND OUTPUT ************************************************************** --> - diff --git a/tests/data/LF_lat_lon_UseCase/prerun_lat_lon.xml b/tests/data/LF_lat_lon_UseCase/prerun_lat_lon.xml index 877fcfbb..29b25053 100644 --- a/tests/data/LF_lat_lon_UseCase/prerun_lat_lon.xml +++ b/tests/data/LF_lat_lon_UseCase/prerun_lat_lon.xml @@ -86,15 +86,6 @@ - - - - diff --git a/tests/data/LF_lat_lon_UseCase/run_lat_lon.xml b/tests/data/LF_lat_lon_UseCase/run_lat_lon.xml index 9ee78767..95765d14 100644 --- a/tests/data/LF_lat_lon_UseCase/run_lat_lon.xml +++ b/tests/data/LF_lat_lon_UseCase/run_lat_lon.xml @@ -87,15 +87,6 @@ - - - - diff --git a/tests/test_caching.py b/tests/test_caching.py index e887c6d5..a2c5c1a7 100644 --- a/tests/test_caching.py +++ b/tests/test_caching.py @@ -16,6 +16,9 @@ class TestCaching(ETRS89TestCase): case_dir = os.path.join(os.path.dirname(__file__), 'data', 'LF_ETRS89_UseCase') settings_file = os.path.join(case_dir, 'settings', 'full.xml') + path_out = os.path.join(case_dir, 'out') + if not os.path.exists(path_out): + os.mkdir(path_out) out_dir_a = os.path.join(case_dir, 'out', 'a') out_dir_b = os.path.join(case_dir, 'out', 'b')