diff --git a/docs/inversion/images/cell_size_estimate.png b/docs/inversion/images/cell_size_estimate.png new file mode 100644 index 00000000..5b104b62 Binary files /dev/null and b/docs/inversion/images/cell_size_estimate.png differ diff --git a/docs/inversion/images/dc_mesh_diagram.png b/docs/inversion/images/dc_mesh_diagram.png new file mode 100644 index 00000000..82676ea9 Binary files /dev/null and b/docs/inversion/images/dc_mesh_diagram.png differ diff --git a/docs/inversion/images/extent_to_pad.png b/docs/inversion/images/extent_to_pad.png new file mode 100644 index 00000000..d8dbdafd Binary files /dev/null and b/docs/inversion/images/extent_to_pad.png differ diff --git a/docs/inversion/images/mag_mesh_diagram.png b/docs/inversion/images/mag_mesh_diagram.png new file mode 100644 index 00000000..170896a1 Binary files /dev/null and b/docs/inversion/images/mag_mesh_diagram.png differ diff --git a/docs/inversion/images/mag_mesh_options.png b/docs/inversion/images/mag_mesh_options.png new file mode 100644 index 00000000..2310e1d7 Binary files /dev/null and b/docs/inversion/images/mag_mesh_options.png differ diff --git a/docs/inversion/images/mesh_core.png b/docs/inversion/images/mesh_core.png index 76501472..eb27bee7 100644 Binary files a/docs/inversion/images/mesh_core.png and b/docs/inversion/images/mesh_core.png differ diff --git a/docs/inversion/images/mesh_mag.png b/docs/inversion/images/mesh_mag.png deleted file mode 100644 index e4a3207e..00000000 Binary files a/docs/inversion/images/mesh_mag.png and /dev/null differ diff --git a/docs/inversion/mesh_design.ipynb b/docs/inversion/mesh_design.ipynb index 23c54a9f..43c1674a 100644 --- a/docs/inversion/mesh_design.ipynb +++ b/docs/inversion/mesh_design.ipynb @@ -9,34 +9,23 @@ "\n", "# Mesh design\n", "\n", - "This section provides details about the mesh creation towards joint inversions. The goal is to create inversion meshes that are well adapted to each survey, while also optimal for the subsequent joint inversion process.\n", - "\n", - "For further details and examples on how to create octree meshes, please visit the [Octree Mesh Creation](https://mirageoscience-octree-creation-app.readthedocs-hosted.com/en/latest) documentation provided by [Mira Geoscience](https://www.mirageoscience.com/mining-industry-software/python-integration).\n", + "An appropriately designed computational mesh is an important component to the geophysical inversion. There are both general and survey-specific strategies that may be used to design 'good' meshes for inversion. This section provides complementary details to the core [Octree creation](https://mirageoscience-octree-creation-app.readthedocs-hosted.com/en/latest) documentation, so readers may want to reference both. \n", "\n", "\n", "## Base parameters\n", "\n", - "In a first step, we need to decide on\n", + "An Octree mesh can be though of as the combination of a base mesh and refinements. The base mesh is defined by the core extent, padding, and the base cell resolution. These concepts are explored in further detail in the following sub-sections.\n", "\n", "- [A region of interest (extent)](#Mesh-extent)\n", "- [Base cell resolution](#Base-cell-size)\n", "- [Padding distance](#Padding-cells)\n", "\n", - "These parameters define our core mesh parameters. While each survey can have its own discretization (refinements), the joint inversion routine requires all meshes to have the same outer extent and base cell size. The following sections provide details on the choices made.\n", "\n", + "![base_mesh](./images/mesh_core.png)\n", "\n", "### Mesh extent\n", "\n", - "Since we need our three simulations (direct-current, magnetics and gravity) to have the same core parameters, we need to create a \"combined\" survey that contains all source and receiver locations in our domain.\n", - "\n", - "\n", - "```{figure} ./images/common_survey.png\n", - "---\n", - "scale: 50%\n", - "name: common_survey\n", - "---\n", - "Creation of a combined curve object from all surveys.\n", - "```\n", + "The mesh extent is provided as a `geoh5py.object` from which the bounding box is computed. For geophysical inversion this will in most cases be the geophysical survey object so that core region of the mesh will be centered below the data. For all `geoh5py.BaseSurvey` objects that contain a `complement` object (transmitters, base stations, current electrodes), the extent will be computed from the superposition of the object and it's complement.\n", "\n", "### Base cell size\n", "\n", @@ -46,16 +35,14 @@ "- Terrain clearance (height above ground)\n", "- Modeling resolution (size of the target)\n", "\n", - "Numerical artefacts due to the discretization of topography decay sufficiently when receivers are roughly two cell dimensions away from the nearest edge or corner. Since the lowest drape height from the airborne magnetics is 60 m, a prudent cell size would be in the range of ~30 m. We also want to avoid computing electrical potential differences within the same cell for numerical accuracy. This is condition is guaranteed for cell dimensions less than half the smallest dipole separation - suggesting an even smaller cell (<=20 m) dimension for the 40m dipole direct-current data simulated here.\n", + "Numerical artefacts due to the discretization of topography decay sufficiently when receivers are roughly two cell dimensions away from the nearest edge or corner of a mesh cell in the air/earth boundary. For example, if the lowest drape height in given airborne survey is 60 m, a prudent cell size would be in the range of ~30 m. Another rule of thumb is that the cell size should be smaller than half the smallest data spacing. This is especially important for electrical methods where the electric field may be discontinuous and numerical artefacts may arise if the cells are not small enough to capture rapid changes in the field. For example, a 40m dipole direct-current survey would suggest a cell size around ~20m.\n", "\n", "\n", "### Padding cells\n", "\n", - "The area of interest covers roughly 2 km in width. As a general rule of thumb, the padding region should be at least as wide as the data extent in order to easily model features with wavelengths that may extend beyond the surveyed area.\n", - "\n", - "(diffusion-distance)=\n", + "As a general rule of thumb, the padding region should be at least as wide as the data extent in order to easily model features with wavelengths that may extend beyond the surveyed area.\n", "\n", - "In the case of EM modeling, we also need to consider the diffusion distance of the EM fields. The [skin depth](http://em.geosci.xyz/content/maxwell1_fundamentals/harmonic_planewaves_homogeneous/skindepth.html?highlight=skin%20depth#approximations) can be estimated by\n", + "In the case of EM modeling, we also need to consider the diffusion distance of the EM fields. The [skin depth](http://em.geosci.xyz/content/maxwell1_fundamentals/harmonic_planewaves_homogeneous/skindepth.html?highlight=skin%20depth#approximations) represents the distance over which the fields decay by a factor of $1/e$. The skin depth can thus be used to add a padding distance that incorporates the minimum area within which the fields have any influence on the solution. The skin depth can be calculated for the frequency-domain system by:\n", "\n", "$$\n", "\\delta = 503 \\sqrt{\\frac{\\rho}{f}}\n", @@ -63,7 +50,7 @@ "\n", "where $\\rho$ is the expected resistivity of the background and $f$ is the base frequency of the system.\n", "\n", - "Equivalently for time-domain systems\n", + "Equivalently, for time-domain systems\n", "\n", "\n", "$$\n", @@ -73,20 +60,7 @@ "where $t\\;(sec)$ is the largest time measured by the system, $\\rho \\; (\\Omega.m)$ is the expected resistivity of the background and $\\mu_0 \\; (4 \\pi * 1e-7)$ is the permeability of free space.\n", "\n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "```{figure} ./images/mesh_core.png\n", - "---\n", - "scale: 50%\n", - "\n", - "name: mesh_core\n", - "---\n", - "Core parameters used for the mesh creation of each inversion. Only the refinement strategies differ between surveys.\n", - "```" + "\n" ] }, { @@ -96,7 +70,7 @@ "source": [ "## Refinements\n", "\n", - "An octree mesh allows to increase the mesh resolution (smaller cells) in specific regions - warranted for either numerical accuracy or modeling purposes. Fine cells can be \"added\" to the octree grid using various insertion methods that depend on the type of object used:\n", + "An octree mesh can be refined in specific regions - warranted for either numerical accuracy or modeling purposes. Fine cells can be \"added\" to the octree grid using various insertion methods that depend on the type of object used:\n", "\n", "- Radial: Add concentric spherical shells of cells around points. It is the default behaviour for `Points`-like objects.\n", "- Line: Add concentric tubes around line segments. It is the default discretization for line sources (EM loops) and `Curve` objects.\n", @@ -111,30 +85,44 @@ "\n", "For EM methods, such as electrical direct-current surveys, the numerical accuracy of the forward simulation is the primary factor to decide on a mesh refinement strategy. It is important to have several small cells around each source and receivers for solving partial differential equations accurately. In this case, dipole source and receivers must be converted to a `Points` object in order to trigger a `radial` refinement. Concentric shells of cells are added around each source/receiver poles.\n", "\n", - "```{figure} ./images/mesh_dc.png\n", - "---\n", - "scale: 50%\n", - "name: mesh_dc\n", - "---\n", - "Parameters used for the mesh creation of the direct-current simulation.\n", - "```\n", + "![dc_mesh](./images/dc_mesh_diagram.png)\n", "\n", - "A second level of `Surface` refinement was added along the triangulated topography. Only coarse cells (3th octree level) were needed to preserve good continuity of the fields near the boundary of the domain.\n", "\n", + "A second level of `Surface` refinement was added along the triangulated topography. Only coarse cells (3th octree level) were needed to preserve good continuity of the fields near the boundary of the domain.\n", "\n", "### For gravity and magnetic surveys\n", "\n", "\n", "For potential field methods, the forward simulation is achieved via integration over prisms (linear operator) which is much less sensitive to the choice of discretization. Moreover, only cells below topography (active cells) are considered. For a more optimal mesh design, only fine cells are required at the air-ground interface. A `horizon` style refinement was used in order to add layers of fine cells below the footprint of the gravity and magnetic sensors. The maximum `distance` parameter was set to 100 m to reduce the number of small cells away from the receiver locations.\n", "\n", - "```{figure} ./images/mesh_mag.png\n", - "---\n", - "scale: 50%\n", - "name: mesh_mag\n", + "![mag_mesh](./images/mag_mesh_diagram.png)" + ] + }, + { + "cell_type": "markdown", + "id": "68f686f7-6458-4b45-b5ff-39c437ba6a97", + "metadata": {}, + "source": [ + "## Auto-mesh\n", + "\n", + "The mesh can be generated in an automatic fashion where the meshing parameters are estimated from the geometry of survey and topography. The mesh parameter in the SimPEG ui.json files is optional. If the user leaves the mesh empty, a mesh will be created during setup for the subsequent inversion. This section describes how the meshing parameters are estimated for the auto-mesh routine.\n", + "\n", + "### Core parameters\n", "\n", - "---\n", - "Parameters used for the mesh creation of the EM simulation.\n", - "```\n" + "The auto-mesh routine refines radially about the survey locations and below the topography surface. The refinement for the survey creates three concentric shells of four cells at each level. The topography refinement places a single layer of two cells in the third octree level below the earth-air interface. The minimum level is set to six adding a base refinement to the whole mesh so that the cells in the padding region do not coarsen excessively.\n", + "\n", + "### Cell size\n", + "\n", + "The cell size will be chosen as half the median station spacing of the survey.\n", + "\n", + "![cell_size](./images/cell_size_estimate.png)\n", + "\n", + "### Padding\n", + "\n", + "The padding size can be estimated from the extent of the survey object. The auto-mesh will first compute the horizontal extent and take half of the largest extent (in x/y) for the horizontal padding parameters.\n", + "\n", + "![extent_to_pad](./images/extent_to_pad.png)\n", + "\n" ] } ], @@ -143,6 +131,18 @@ "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.16" } }, "nbformat": 4, diff --git a/docs/inversion/mesh_design.py b/docs/inversion/mesh_design.py index 8508fe05..319ef760 100644 --- a/docs/inversion/mesh_design.py +++ b/docs/inversion/mesh_design.py @@ -16,34 +16,23 @@ # # # Mesh design # -# This section provides details about the mesh creation towards joint inversions. The goal is to create inversion meshes that are well adapted to each survey, while also optimal for the subsequent joint inversion process. -# -# For further details and examples on how to create octree meshes, please visit the [Octree Mesh Creation](https://mirageoscience-octree-creation-app.readthedocs-hosted.com/en/latest) documentation provided by [Mira Geoscience](https://www.mirageoscience.com/mining-industry-software/python-integration). +# An appropriately designed computational mesh is an important component to the geophysical inversion. There are both general and survey-specific strategies that may be used to design 'good' meshes for inversion. This section provides complementary details to the core [Octree creation](https://mirageoscience-octree-creation-app.readthedocs-hosted.com/en/latest) documentation, so readers may want to reference both. # # # ## Base parameters # -# In a first step, we need to decide on +# An Octree mesh can be though of as the combination of a base mesh and refinements. The base mesh is defined by the core extent, padding, and the base cell resolution. These concepts are explored in further detail in the following sub-sections. # # - [A region of interest (extent)](#Mesh-extent) # - [Base cell resolution](#Base-cell-size) # - [Padding distance](#Padding-cells) # -# These parameters define our core mesh parameters. While each survey can have its own discretization (refinements), the joint inversion routine requires all meshes to have the same outer extent and base cell size. The following sections provide details on the choices made. # +# ![base_mesh](./images/mesh_core.png) # # ### Mesh extent # -# Since we need our three simulations (direct-current, magnetics and gravity) to have the same core parameters, we need to create a "combined" survey that contains all source and receiver locations in our domain. -# -# -# ```{figure} ./images/common_survey.png -# --- -# scale: 50% -# name: common_survey -# --- -# Creation of a combined curve object from all surveys. -# ``` +# The mesh extent is provided as a `geoh5py.object` from which the bounding box is computed. For geophysical inversion this will in most cases be the geophysical survey object so that core region of the mesh will be centered below the data. For all `geoh5py.BaseSurvey` objects that contain a `complement` object (transmitters, base stations, current electrodes), the extent will be computed from the superposition of the object and it's complement. # # ### Base cell size # @@ -53,16 +42,14 @@ # - Terrain clearance (height above ground) # - Modeling resolution (size of the target) # -# Numerical artefacts due to the discretization of topography decay sufficiently when receivers are roughly two cell dimensions away from the nearest edge or corner. Since the lowest drape height from the airborne magnetics is 60 m, a prudent cell size would be in the range of ~30 m. We also want to avoid computing electrical potential differences within the same cell for numerical accuracy. This is condition is guaranteed for cell dimensions less than half the smallest dipole separation - suggesting an even smaller cell (<=20 m) dimension for the 40m dipole direct-current data simulated here. +# Numerical artefacts due to the discretization of topography decay sufficiently when receivers are roughly two cell dimensions away from the nearest edge or corner of a mesh cell in the air/earth boundary. For example, if the lowest drape height in given airborne survey is 60 m, a prudent cell size would be in the range of ~30 m. Another rule of thumb is that the cell size should be smaller than half the smallest data spacing. This is especially important for electrical methods where the electric field may be discontinuous and numerical artefacts may arise if the cells are not small enough to capture rapid changes in the field. For example, a 40m dipole direct-current survey would suggest a cell size around ~20m. # # # ### Padding cells # -# The area of interest covers roughly 2 km in width. As a general rule of thumb, the padding region should be at least as wide as the data extent in order to easily model features with wavelengths that may extend beyond the surveyed area. -# -# (diffusion-distance)= +# As a general rule of thumb, the padding region should be at least as wide as the data extent in order to easily model features with wavelengths that may extend beyond the surveyed area. # -# In the case of EM modeling, we also need to consider the diffusion distance of the EM fields. The [skin depth](http://em.geosci.xyz/content/maxwell1_fundamentals/harmonic_planewaves_homogeneous/skindepth.html?highlight=skin%20depth#approximations) can be estimated by +# In the case of EM modeling, we also need to consider the diffusion distance of the EM fields. The [skin depth](http://em.geosci.xyz/content/maxwell1_fundamentals/harmonic_planewaves_homogeneous/skindepth.html?highlight=skin%20depth#approximations) represents the distance over which the fields decay by a factor of $1/e$. The skin depth can thus be used to add a padding distance that incorporates the minimum area within which the fields have any influence on the solution. The skin depth can be calculated for the frequency-domain system by: # # $$ # \delta = 503 \sqrt{\frac{\rho}{f}} @@ -70,7 +57,7 @@ # # where $\rho$ is the expected resistivity of the background and $f$ is the base frequency of the system. # -# Equivalently for time-domain systems +# Equivalently, for time-domain systems # # # $$ @@ -82,22 +69,10 @@ # # # -# -# -# -# -# ```{figure} ./images/mesh_core.png -# --- -# scale: 50% -# -# name: mesh_core -# --- -# Core parameters used for the mesh creation of each inversion. Only the refinement strategies differ between surveys. -# ``` # ## Refinements # -# An octree mesh allows to increase the mesh resolution (smaller cells) in specific regions - warranted for either numerical accuracy or modeling purposes. Fine cells can be "added" to the octree grid using various insertion methods that depend on the type of object used: +# An octree mesh can be refined in specific regions - warranted for either numerical accuracy or modeling purposes. Fine cells can be "added" to the octree grid using various insertion methods that depend on the type of object used: # # - Radial: Add concentric spherical shells of cells around points. It is the default behaviour for `Points`-like objects. # - Line: Add concentric tubes around line segments. It is the default discretization for line sources (EM loops) and `Curve` objects. @@ -112,28 +87,36 @@ # # For EM methods, such as electrical direct-current surveys, the numerical accuracy of the forward simulation is the primary factor to decide on a mesh refinement strategy. It is important to have several small cells around each source and receivers for solving partial differential equations accurately. In this case, dipole source and receivers must be converted to a `Points` object in order to trigger a `radial` refinement. Concentric shells of cells are added around each source/receiver poles. # -# ```{figure} ./images/mesh_dc.png -# --- -# scale: 50% -# name: mesh_dc -# --- -# Parameters used for the mesh creation of the direct-current simulation. -# ``` +# ![dc_mesh](./images/dc_mesh_diagram.png) # -# A second level of `Surface` refinement was added along the triangulated topography. Only coarse cells (3th octree level) were needed to preserve good continuity of the fields near the boundary of the domain. # +# A second level of `Surface` refinement was added along the triangulated topography. Only coarse cells (3th octree level) were needed to preserve good continuity of the fields near the boundary of the domain. # # ### For gravity and magnetic surveys # # # For potential field methods, the forward simulation is achieved via integration over prisms (linear operator) which is much less sensitive to the choice of discretization. Moreover, only cells below topography (active cells) are considered. For a more optimal mesh design, only fine cells are required at the air-ground interface. A `horizon` style refinement was used in order to add layers of fine cells below the footprint of the gravity and magnetic sensors. The maximum `distance` parameter was set to 100 m to reduce the number of small cells away from the receiver locations. # -# ```{figure} ./images/mesh_mag.png -# --- -# scale: 50% -# name: mesh_mag +# ![mag_mesh](./images/mag_mesh_diagram.png) + +# ## Auto-mesh +# +# The mesh can be generated in an automatic fashion where the meshing parameters are estimated from the geometry of survey and topography. The mesh parameter in the SimPEG ui.json files is optional. If the user leaves the mesh empty, a mesh will be created during setup for the subsequent inversion. This section describes how the meshing parameters are estimated for the auto-mesh routine. +# +# ### Core parameters +# +# The auto-mesh routine refines radially about the survey locations and below the topography surface. The refinement for the survey creates three concentric shells of four cells at each level. The topography refinement places a single layer of two cells in the third octree level below the earth-air interface. The minimum level is set to six adding a base refinement to the whole mesh so that the cells in the padding region do not coarsen excessively. +# +# ### Cell size +# +# The cell size will be chosen as half the median station spacing of the survey. +# +# ![cell_size](./images/cell_size_estimate.png) +# +# ### Padding +# +# The padding size can be estimated from the extent of the survey object. The auto-mesh will first compute the horizontal extent and take half of the largest extent (in x/y) for the horizontal padding parameters. +# +# ![extent_to_pad](./images/extent_to_pad.png) # -# --- -# Parameters used for the mesh creation of the EM simulation. -# ``` #