Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Create bathymetry mesh + refinement options #301

Closed
MatiasAlday opened this issue Feb 22, 2022 · 16 comments
Closed

Create bathymetry mesh + refinement options #301

MatiasAlday opened this issue Feb 22, 2022 · 16 comments

Comments

@MatiasAlday
Copy link

MatiasAlday commented Feb 22, 2022

Hello everyone,

I'm new to Thetis, so far I've run the examples and played a bit with them. I have many questions but I will leave 4 for the moment:
First the context: I need to set up a tidal model

  1. How can I create an unstructured (triangles) mesh which domain is not rectangular (as in the examples) and where I need to apply some refinement areas, most likely close to the shore. ? Or have an automatic triangle size change based on depth and/or CFL ?

  2. How should I define the coastline and the open boundaries ? (I'm using OpenStreetmap polygons and Emodnet Bathymetry)

  3. After creating the mesh (if I succeed), how can I access the coordinates of each node of the bathymetry mesh and the corresponding nodes defining an element (triangle), etc. ? This can be generalized for all field output since I want to be able to plot and handle the data with my own post process script outside Paraview.

  4. any recommendations on how to handle data after running a model? I understand that it's not advice to mess with the firedrake environment. So I was thinking maybe loading the default output .pvd files and working with them in a different environment.

@jhill1
Copy link
Contributor

jhill1 commented Feb 22, 2022 via email

@MatiasAlday
Copy link
Author

Hi Jon ,
Thanks a lot for the input and the quick reply.
I'll start with checking gmesh and generating something there. Then, as far as I understand, I should actually use qmesh to create the initial version of the mesh with coastlines and open boundaries defined, and then apply refinement in gmesh.

cheers,
Matias

@jwallwork23
Copy link
Contributor

Regarding the mesh coordinates, they are represented as a Firedrake Function in the mesh's coordinates attribute. In the special (but fairly common in Thetis codes) case that both your mesh and the bathymetry field are P1, mesh.coordinates.dat.data gives you the vertex array on the current process and bathymetry.dat.data gives the corresponding values. Is this the case, or are you using other finite elements?

@MatiasAlday
Copy link
Author

MatiasAlday commented Feb 23, 2022

Yes, that I could realize. But for example the nodes (or vertex) set defining a mesh element I could not find. Or in the particular case of the 3rd tutorial about a 3D channel with boundary conditions, I couldn't find the field values related to each layer.

@jwallwork23
Copy link
Contributor

jwallwork23 commented Feb 23, 2022

One way to do this is to use the PETSc DMPlex mesh data structure that underpins Firedrake's mesh concept. The plex object has a different numbering to the Firedrake mesh and uses unique indices across all entities.

First, make the numberings consistent between the two mesh concepts as described here.

Then something like this should do the trick in the non-extruded case:

plex = mesh.topology_dm

vStart, vEnd = plex.getDepthStratum(0)
vertices = set(range(vStart, vEnd))  # set of all vertex indices

cStart, cEnd = plex.getHeightStratum(0)
for c in range(cStart, cEnd):  # loop over all cell indices
    closure = set(plex.getTransitiveClosure(c)[0])  # get all entities associated with the cell
    vloc = closure.intersection(vertices)  # restrict to vertices associated with the cell (disregard edges)
    print(f"Cell {c} has vertex set {vloc}")

The doc page I linked above tells you how to get the coordinate values for those vertices.

All of this is probably also possible using a par_loop.

Hopefully this is useful!

@MatiasAlday
Copy link
Author

Thanks! I'll take a look at it . I'll let you know if I have something (or can't get anything ).
Cheers,
Matias

@tkarna
Copy link
Contributor

tkarna commented Feb 23, 2022

In order to evaluate model fields at points, I recommend storing time series at run time:

tscb = TimeSeriesCallback2D(['elev_2d'], x, y, 'StationA', append_to_log=False)
solver_obj.add_callback(tscb)

where x and y are in mesh coordinates. This will create a hdf5 file outputs/diagnostic_timeseries_StationA_elev.hdf5. The time is model time in seconds.

In case of a 3D model, you can also extract vertical profiles (VerticalProfileCallback) and transects (TransectCallback).

@MatiasAlday
Copy link
Author

Ah great, so then I can just load and handle the stored data in the hdf5 file. I'll come back to you once I learn all of the bathy mesh generation up to reading the mesh coordinates part.

@tkarna
Copy link
Contributor

tkarna commented Feb 23, 2022

At the moment we do not have portable output files that would allow easy post processing (for experimental netcdf outputs see #255).

If you wish to plot 2D fields in post-processing this can be done in Python by reading in the vtk output files, e.g.

import matplotlib
import matplotlib.pyplot as plt
import vtk
from vtk.util.numpy_support import vtk_to_numpy

def read_vtu(f, parrayname='Salinity', values_only=False, verbose=True,
             parallel=False):
    """
    Read unstructured 2D mesh from pvtu file to numpy arrays.

    Assumes a triangular mesh.
    """
    if verbose:
        print('Reading {:}'.format(f))
    if parallel:
        reader = vtk.vtkXMLPUnstructuredGridReader()
    else:
        reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(f)
    reader.Update()

    data = reader.GetOutput()

    xyz = elems = None
    if not values_only:
        # read grid
        points = data.GetPoints()
        xyz = vtk_to_numpy(points.GetData())
        elems = vtk_to_numpy(data.GetCells().GetData())
        # NOTE this depends on the element type and polynomial degree
        elems = elems.reshape((-1, 4))[:, 1:]

    assert reader.GetPointArrayName(0) == parrayname
    values = vtk_to_numpy(data.GetPointData().GetArray(0))
    return xyz, elems, values
    
# pvtu for parallel outputs
f = 'outputs/Elevation2d/Elevation2d_130.pvtu'
# you can check the array name in the vtu/pvtu file
xyz, elems, values = read_vtu(f, parrayname='Elevation', parallel=True)
tri = matplotlib.tri.Triangulation(xyz[:, 0], xyz[:, 1], elems)
plt.tripcolor(tri, values, ...)

@MatiasAlday
Copy link
Author

By the way , is it OK if I pip install matplotlib and other libraries in the firedrake env? or is it better to have separate env to play with the output ?

@tkarna
Copy link
Contributor

tkarna commented Feb 23, 2022

is it OK if I pip install matplotlib and other libraries in the firedrake env?

That's fine. I think only anaconda is incompatible, so in that case it's better to have separate python environments

@MatiasAlday
Copy link
Author

Hi Matias, Replies inline below.
On 22/02/2022 09:52, MatiasAlday wrote: Hello everyone, I'm new to Thetis, so far I've run the examples and played a bit with them. I have many questions but I will leave 4 for the moment: First the context: I need to set up a tidal model 1. How can I create an unstructured (triangles) mesh which domain is not rectangular (as in the examples) and where I need to apply some refinement areas, most likely close to the shore. ? Or have an automatic triangle size change based on depth and/or CFL ?
You can use the refinement via gmsh. You can create any metric you wish (e.g. CFL or depth) via loading in an FLD file as a structured grid and then using gmsh's maths routines. An example .geo file is attached for gmsh v2-something. FLDs can be created using qmesh (see below)
2. How should I define the coastline and the open boundaries ? (I'm using OpenStreetmap polygons and Emodnet Bathymetry)
I would recommend qmesh. This is currently in development for QGIS3 python3, but does work when installed manually using python3. https://bitbucket.org/qmesh-developers/workspace/projects/QMES You'll need to install the setuptools-qmesh, then the main qmesh, then the qmesh-cli. I currently use this to make all meshes. Note I'm an author on the package and others are available!
3. After creating the mesh (if I succeed), how can I access the coordinates of each node of the bathymetry mesh and the corresponding nodes defining an element (triangle), etc. ? This can be generalized for al field output since I want to be able to plot and handle the data with my own pot process script outside Paraview.
I'll leave others to answer this in more detail, but you can access this via thetis/numpy type code.
4. any recommendations on how to handle data after running a model? I understand that it's not advice to mess with the firedrake environment. So I was thinking maybe loading the default output .pvd files and working with them in a different environment.
I use paraview and/or python to extract the data I need. Hope that helps, Jon

By the way, does qmesh work only on Linux based OS? or can it be installed on the Windows version of python?

@jhill1
Copy link
Contributor

jhill1 commented Feb 23, 2022 via email

@MatiasAlday
Copy link
Author

MatiasAlday commented Feb 23, 2022

Jon, sorry to annoy you with basic questions.
I keep getting this message -> ERROR: Failed building wheel for qmesh
The same happens with GDAL .

I tried with pip install qmesh --no-cache-dir just in case, but I keep getting the error.
Any Ideas?

Added: Actually I can't build the wheel for the latest GDAL version . When I execute pip install qmesh , it points me to the latest version of GDAL. I just upgraded to v3.3.3 and worked just fine:
(mesh_env) maldaygonzalez@TUD260721:~$ pip install GDAL==3.3.3
Collecting GDAL==3.3.3
Downloading GDAL-3.3.3.tar.gz (747 kB)
|████████████████████████████████| 747 kB 2.1 MB/s
Building wheels for collected packages: GDAL
Building wheel for GDAL (setup.py) ... done
Created wheel for GDAL: filename=GDAL-3.3.3-cp38-cp38-linux_x86_64.whl size=3158054 sha256=01a7132016ee698b5fe94a53ce078f7bfdf789e4eeedefe00273fe69c04f442f
Stored in directory: /home/maldaygonzalez/.cache/pip/wheels/67/de/db/3e3c736860be40b33ca94c9c727731ae09294c20d815cafbc3
Successfully built GDAL
Installing collected packages: GDAL
Successfully installed GDAL-3.3.3

@jhill1
Copy link
Contributor

jhill1 commented Feb 24, 2022 via email

@AlexandrosAvdis
Copy link

AlexandrosAvdis commented Mar 1, 2022 via email

@tkarna tkarna closed this as completed May 20, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants