grid.rst 6.34 KB
Newer Older
1 2 3
Grid Creation and Mapping
=========================

Lukas Riedel's avatar
Lukas Riedel committed
4 5
To guarantee numeric accuracy, :doc:`boundary conditions <bcfile>`
and :doc:`parameterizations <parameter-file>` must exactly map to grid
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
boundary faces and grid cells, respectively. DORiE therefore not only requires
a specification of the actual grid, but also of these mappings.

DORiE supports two types of grid input, which are controlled via the config
file parameter ``grid.gridType``:

1. Building a rectangular grid on the fly (``rectangular``).
   Depending on whether grid adaptivity is enabled or not, DORiE will use
   a structured or an unstructured grid manager. This does not change the input
   scheme. The user has to specify the config parameters for ``[grid]``
   ``dimensions``, spatial ``extensions``, and grid ``cells`` in each
   direction. For the mapping of boundary conditions and parameters,
   :ref:`mapping datasets <man-grid_mapping-dataset>` are required.

2. Building a grid according to a GMSH grid file (``gmsh``).
   This mode always uses an unstructured grid manager. The user only has to
   specify the ``[grid]`` ``dimensions`` and the GMSH ``gridFile``.
   The mapping is read from this file and must be considered when
   :ref:`building the mesh <man-grid_gmsh>`.


.. _man-grid_mapping-dataset:

Mapping Datasets
----------------
Mappings are realized as HDF5_ datasets, where every value is mapped to a
single grid entity. Changing the grid configuration requires adapting these
datasets! All config keys of this part refer to the config section
``[grid.mapping]``. There is one dataset for the mapping of the soil
architecture (``volume``) and one dataset for each segment of the boundary
(``boundaryXYZ``). The datasets contain *integer* indices, which are then
identified with certain boundary conditions or parameterizations given by
the respective YAML input files. All mapping datasets must be stored inside
the same H5 file (``file``). Use the Python module ``h5py`` for easily
creating these data files.

The dataset layout follows the C memory layout, with the primary dimension
running fastest. The coordinate system is (right-hand) cartesian, with the
x-axis running to the right in both 2D and 3D.
All datasets follow the general layout ``array[z, y, x]``, where unused
dimensions are simply omitted. The first dimension is always the vertical axis.
The ``volume`` dataset must have the same dimension as the grid itself
(``grid.dimensions``), whereas the ``boundaryXYZ`` datasets are interpreted as
projections of the boundary in their respective normal direction, reducing the
dimension by one. The ``boundaryLower`` dataset would have the layout
``array[x]`` for a 2D grid, and ``array[y, x]`` for a 3D grid. Therefore,
``boundaryLower`` and ``boundaryUpper`` are "seen from above"
(3D-layout: ``array[y, x]``), and ``boundaryLeft`` and ``boundaryRight`` are
"seen from the right" (3D-layout: ``array[z, y]``).

The following Python code creates a volume mapping dataset ``layered`` in a
file ``mapping.h5``, where the lower half of the domain maps to medium 0,
and the upper half to medium 1. The 3D domain contains 10 cells in each
direction.

.. code-block:: python

    import numpy as np
    import h5py

    # Create dataset. Mind the data type!
    size = 10
    layered = np.zeros((size, size, size), dtype=np.int_)
    layered[5:, ...] = 1

    # Write dataset to file in 'append' mode
    with h5py.File("mapping.h5", 'a') as file:
        file.create_dataset("layered", data=layered)

Homogeneous Mappings
^^^^^^^^^^^^^^^^^^^^
If the entire volume, or a complete boundary, should be mapped to a single
index, no dataset is required. Instead, you can set the value for the
respective config file parameter to the desired index. If no dataset should be
read at all, set ``file`` to ``none``.

Even with ``file`` set to a valid H5 file, DORiE will *always* try to
interpret the input values for ``grid.mapping`` as integers. If this succeeds,
the value is interpreted as "global" mapping index for the respective part
of the grid. Therefore, **do not use dataset names starting with digits!**

.. _man-grid_gmsh:

Mapping Soil Layers with GMSH
-----------------------------
GMSH_ supports mapping lines, surfaces, and volumes to *physical* entities.
These entities may combine multiple of the aforementioned *geometric*
entities. Physical entities are assigned a non-negative index upon creation.
These indices can be set by the user in any of the GMSH frontends.
Indices may not occur multiple times, even if they are assigned to different
types of physical entities.

The following ``.geo`` GMSH input file creates a rectangular domain that is
split in half. The lower part is mapped to index 1, the upper part to index 2.
Additionally, a different index for every boundary orientation is set. Notice
that the left and right boundaries consist of two *geometric* lines each.
This code can directly be transferred to a Python script using
the pygmsh_ module. It writes a ``.geo`` file while checking for a correct
syntax within your script. It is readily installed into the
Lukas Riedel's avatar
Lukas Riedel committed
105
:doc:`virtual environment <cli>`.
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

.. code-block:: default

    // define geometric entities
    Point(1) = {0, 0, 0, 0.1};
    Point(2) = {2, 0, 0, 0.1};
    Point(3) = {2, 2, 0, 0.1};
    Point(4) = {0, 2, 0, 0.1};
    Point(5) = {0, 1, 0, 0.1};
    Point(6) = {2, 1, 0, 0.1};

    Line(1) = {1, 2};
    Line(2) = {2, 6};
    Line(3) = {6, 3};
    Line(4) = {3, 4};
    Line(5) = {4, 5};
    Line(6) = {5, 1};
    Line(7) = {5, 6};

    Line Loop(1) = {1, 2, -7, 6};
    Plane Surface(1) = {1};
    Line Loop(2) = {7, 3, 4, 5};
    Plane Surface(2) = {2};

    // define physical entities, index in round brackets
    Physical Surface(1) = {1}; // lower
    Physical Surface(2) = {2}; // upper

    // entire set of physical entities must always be defined!
    Physical Line(3) = {1}; // bottom
    Physical Line(4) = {2, 3}; // right
    Physical Line(5) = {4}; // top
    Physical Line(6) = {5, 6}; // left

A ``.geo`` file is the basis for creating the actual mesh in GMSH. You can
load it into the GMSH GUI, or perform the meshing directly using the
`GMSH command line interface
<http://gmsh.info/doc/texinfo/gmsh.html#Non_002dinteractive-mode>`_:

    gmsh <geo-file> -<dim>

Replace ``<geo-file>`` with the appropriate file, and ``dim`` with the
spatial dimension of the intended mesh.

.. _HDF5: https://www.h5py.org/
.. _GMSH: http://gmsh.info/
.. _pygmsh: https://pypi.org/project/pygmsh/