GeoscienceAustralia / anuga_core

AnuGA for the simulation of the shallow water equation
https://anuga.anu.edu.au
Other
193 stars 94 forks source link

How to setup river level in anuga #270

Open chooron opened 1 year ago

chooron commented 1 year ago

I have prepared the elevation data for the basin and the river polygon. However, when I run this code:

domain.set_quantity('stage', expression='elevation + %f' % 3, polygon=river_polygon)

I encounter the following error:


Traceback (most recent call last):
  File "D:\miniconda3\envs\anuga_env\lib\site-packages\IPython\core\interactiveshell.py", line 3508, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-4-3914e5c91e47>", line 1, in <module>
    domain.set_quantity('stage', expression='elevation + %f' % 3, polygon=river_polygon)
  File "D:\miniconda3\envs\anuga_env\lib\site-packages\anuga\shallow_water\shallow_water_domain.py", line 1058, in set_quantity
    Generic_Domain.set_quantity(self, name, *args, **kwargs)
  File "D:\miniconda3\envs\anuga_env\lib\site-packages\anuga\abstract_2d_finite_volumes\generic_domain.py", line 722, in set_quantity
    self.quantities[name].set_values(*args, **kwargs)
  File "D:\miniconda3\envs\anuga_env\lib\site-packages\anuga\abstract_2d_finite_volumes\quantity.py", line 839, in set_values
    raise Exception(msg)
Exception: With a polygon selected, 'set_quantity' must provide the 'numeric' keyword, and it must be a constant value.

Therefore, I want to know how to initialize the water level of the river using a polygon. When I run the code:

domain.set_quantity('stage', 1900, polygon=river_polygon)

It seems to run without errors, but the river stages are definitely inconsistent."

stoiver commented 1 year ago

@chooron, strange that we don't have set_quantity with polygon implemented properly. Must add that to our todo list.

With your simple example try:

domain.set_quantity('stage', numeric=1900, polygon=river_polygon)

I will try to provide a work around for your problem. Should be easy to code.

stoiver commented 1 year ago

@chooron The following should work, though a bit low level.

region = anuga.Region(domain, polygon=river_polygon, expand_polygon=True)
tids = region.indices # ids of triangles intersecting region defined by river polygon

stage_c = domain.quantities['stage'].centroid_values
elev_c  = domain.quantities['elevation'].centroid_values

stage_c[tids] = elev_c[tids] + 3.0
chooron commented 1 year ago

@chooron The following should work, though a bit low level.

region = anuga.Region(domain, polygon=river_polygon, expand_polygon=True)
tids = region.indices # ids of triangles intersecting region defined by river polygon

stage_c = domain.quantities['stage'].centroid_values
elev_c  = domain.quantities['elevation'].centroid_values

stage_c[tids] = elev_c[tids] + 3.0

Thanks ! I have solved the problem by preparing the stage file by using Arcgis, and run the code:

stage = domain.quantities['stage']
stage.set_values(filename=join(file_path, 'level3m.asc'))
chooron commented 1 year ago

Hi @stoiver , I'm not sure if the watershed is too large, the inflow is too small, etc. My flood evolution is not very effective, and the river depths no longer change significantly after one time period of calculation.

stoiver commented 1 year ago

@chooron what do you mean by "one time period of calculation"?

Setting up the mesh, the elevation and stage can take a bit of work. The interaction of the elevation and the mesh can lead to some interesting artefacts, in particular artificial numerical dams.

One trick is to add a large amount of water to the domain (via rate_operator) to fill up all these dams and change the centroid elevation to match. This can be unsatisfying as the elevation does not match your original DEM.

Refining the mesh should help, but you might end up with an unnecessarily large number of triangles and so a slow simulation.

You can adapt your mesh to refine along rivers so that triangulation follows the river as opposed to crossing the river. If you have a look at the anuga-community/anuga_core issues list, in particular https://github.com/anuga-community/anuga_core/issues/12 then there are some links to producing good meshes, which are finer along the river courses. You may still need to have a burn-in procedure of rain over the catchment to "fill" up the rivers