diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 94c49978..15518c52 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -46,6 +46,7 @@ jobs: - geos-geomechanics - geos-mesh - geos-posp + - geos-processing - geos-timehistory - geos-trame - geos-xml-tools @@ -57,8 +58,10 @@ jobs: dependencies: "geos-utils" - package-name: geos-mesh dependencies: "geos-utils geos-geomechanics" - - package-name: geos-posp + - package-name: geos-processing dependencies: "geos-utils geos-mesh geos-geomechanics" + - package-name: geos-posp + dependencies: "geos-utils geos-mesh geos-geomechanics geos-processing" - package-name: pygeos-tools dependencies: "geos-utils geos-mesh" - package-name: geos-timehistory diff --git a/.github/workflows/typing-check.yml b/.github/workflows/typing-check.yml index 4b53634b..0a00276d 100644 --- a/.github/workflows/typing-check.yml +++ b/.github/workflows/typing-check.yml @@ -1,5 +1,5 @@ name: "Pull Request Typing Check" -on: +on: - pull_request # Cancels in-progress workflows for a PR when updated @@ -16,7 +16,7 @@ jobs: max-parallel: 3 matrix: # add packages to check typing - package-name: ["geos-geomechanics", "geos-posp", "geos-timehistory", "geos-utils", "geos-trame", "geos-xml-tools", "hdf5-wrapper"] + package-name: ["geos-geomechanics", "geos-posp", "geos-processing", "geos-timehistory", "geos-utils", "geos-trame", "geos-xml-tools", "hdf5-wrapper"] steps: - uses: actions/checkout@v4 diff --git a/docs/conf.py b/docs/conf.py index 7d925de4..7f9197b0 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -17,7 +17,7 @@ # Add python modules to be documented python_root = '..' -python_modules = ( 'geos-ats', 'geos-geomechanics', 'geos-mesh', 'geos-posp', 'geos-pv', 'geos-timehistory', +python_modules = ( 'geos-ats', 'geos-geomechanics', 'geos-mesh', 'geos-posp', 'geos-processing', 'geos-pv', 'geos-timehistory', 'geos-utils', 'geos-xml-tools', 'geos-xml-viewer', 'hdf5-wrapper', 'pygeos-tools' ) @@ -45,7 +45,7 @@ # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. extensions = [ - 'sphinx.ext.napoleon', 'sphinx_design', 'sphinx_rtd_theme', 'sphinxarg.ext', 'sphinxcontrib.programoutput', + 'sphinx.ext.napoleon', 'sphinx_design', 'sphinx_rtd_theme', 'sphinxarg.ext', 'sphinxcontrib.programoutput', 'sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.mathjax', 'sphinx.ext.todo', 'sphinx.ext.viewcode' ] diff --git a/docs/geos-processing.rst b/docs/geos-processing.rst new file mode 100644 index 00000000..735071d9 --- /dev/null +++ b/docs/geos-processing.rst @@ -0,0 +1,13 @@ +GEOS Processing tools +======================== + + +.. toctree:: + :maxdepth: 5 + :caption: Contents: + + ./geos_processing_docs/post_processing.rst + + ./geos_processing_docs/pre_processing.rst + + ./geos_processing_docs/generic_processing_tools.rst diff --git a/docs/geos_mesh_docs/processing.rst b/docs/geos_mesh_docs/processing.rst index fc8f3ca2..b5037228 100644 --- a/docs/geos_mesh_docs/processing.rst +++ b/docs/geos_mesh_docs/processing.rst @@ -27,7 +27,6 @@ geos.mesh.processing.FillPartialArrays filter :undoc-members: :show-inheritance: - geos.mesh.processing.meshQualityMetricHelpers module ----------------------------------------------------- diff --git a/docs/geos_posp_docs/PVplugins.rst b/docs/geos_posp_docs/PVplugins.rst index f95f450c..8f41728b 100644 --- a/docs/geos_posp_docs/PVplugins.rst +++ b/docs/geos_posp_docs/PVplugins.rst @@ -34,13 +34,6 @@ PVExtractMergeBlocksVolumeWell plugin .. automodule:: PVplugins.PVExtractMergeBlocksVolumeWell - -PVGeomechanicsAnalysis plugin ---------------------------------------- - -.. automodule:: PVplugins.PVGeomechanicsAnalysis - - PVGeomechanicsWorkflowVolume plugin --------------------------------------------- diff --git a/docs/geos_posp_docs/filters.rst b/docs/geos_posp_docs/filters.rst index 30ec1c8d..b3c3cd89 100644 --- a/docs/geos_posp_docs/filters.rst +++ b/docs/geos_posp_docs/filters.rst @@ -3,14 +3,6 @@ vtk Filters This package defines vtk filters that allows to process Geos outputs. -geos_posp.filters.GeomechanicsCalculator module ---------------------------------------------------- - -.. automodule:: geos_posp.filters.GeomechanicsCalculator - :members: - :undoc-members: - :show-inheritance: - geos_posp.filters.GeosBlockExtractor module ----------------------------------------------- diff --git a/docs/geos_processing_docs/generic_processing_tools.rst b/docs/geos_processing_docs/generic_processing_tools.rst new file mode 100644 index 00000000..4fdb32fc --- /dev/null +++ b/docs/geos_processing_docs/generic_processing_tools.rst @@ -0,0 +1,4 @@ +Generic processing filters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In progress.. \ No newline at end of file diff --git a/docs/geos_processing_docs/post_processing.rst b/docs/geos_processing_docs/post_processing.rst new file mode 100644 index 00000000..b376bdbf --- /dev/null +++ b/docs/geos_processing_docs/post_processing.rst @@ -0,0 +1,15 @@ +Post-processing filters +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This module contains GEOS post-processing tools. + +Geomechanics post-processing tools +===================================== + +geos.processing.post_processing.GeomechanicsCalculator module +-------------------------------------------------------------- + +.. automodule:: geos.processing.post_processing.GeomechanicsCalculator + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/geos_processing_docs/pre_processing.rst b/docs/geos_processing_docs/pre_processing.rst new file mode 100644 index 00000000..85973093 --- /dev/null +++ b/docs/geos_processing_docs/pre_processing.rst @@ -0,0 +1,4 @@ +Pre-processing filters +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In progress.. \ No newline at end of file diff --git a/docs/geos_pv_docs/processing.rst b/docs/geos_pv_docs/processing.rst index 42fb5590..7b3753b6 100644 --- a/docs/geos_pv_docs/processing.rst +++ b/docs/geos_pv_docs/processing.rst @@ -6,6 +6,7 @@ PVAttributeMapping .. automodule:: geos.pv.plugins.PVAttributeMapping + PVCreateConstantAttributePerRegion ----------------------------------- @@ -18,6 +19,12 @@ PVFillPartialArrays .. automodule:: geos.pv.plugins.PVFillPartialArrays +PVGeomechanicsCalculator plugin +--------------------------------------- + +.. automodule:: geos.pv.plugins.PVGeomechanicsCalculator + + PVSplitMesh ---------------------------------- diff --git a/docs/index.rst b/docs/index.rst index fedfdeb2..b7175f8d 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -47,7 +47,7 @@ To do this, you can clone a copy of the geosPythonPackages repository and instal cd geosPythonPackages/ python -m pip install --upgrade geos-ats - + .. note:: To upgrade an existing installation, the python executable in the above command should correspond to the version you indicated in your host config. If you have previously built the tools, this version will be linked in the build directory: `build_dir/bin/python`. @@ -87,7 +87,9 @@ Packages geos-mesh geos-posp - + + geos-processing + geos-pv geos-timehistory diff --git a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py index 44aa940b..4c7b16e7 100644 --- a/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py +++ b/geos-geomechanics/src/geos/geomechanics/processing/geomechanicsCalculatorFunctions.py @@ -506,8 +506,7 @@ def reservoirStressPathReal( deltaStress: npt.NDArray[ np.float64 ], den: npt.NDArray[ np.float64 ] = np.copy( deltaPressure ) den[ mask ] = 1.0 # use -1 to agree with Geos convention (i.e., compression with negative stress) - # take the xx, yy, and zz components only - rsp: npt.NDArray[ np.float64 ] = np.copy( deltaStress[ :, :3 ] ) + rsp: npt.NDArray[ np.float64 ] = np.copy( deltaStress[ :, : ] ) for j in range( rsp.shape[ 1 ] ): rsp[ :, j ] /= den rsp[ mask, j ] = np.nan diff --git a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py index 8ea4f11d..7fc8b4b5 100644 --- a/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py +++ b/geos-geomechanics/tests/testsFunctionsGeomechanicsCalculator.py @@ -181,11 +181,11 @@ def test_ReservoirStressPathReal( self: Self ) -> None: deltaPressure: npt.NDArray[ np.float64 ] = 0.2 * pressure obtained: npt.NDArray[ np.float64 ] = fcts.reservoirStressPathReal( deltaStress, deltaPressure ) expected: npt.NDArray[ np.float64 ] = np.array( [ - [ 8.0, 12.5, 9.0 ], - [ 0.32, 0.46, 0.34 ], - [ np.nan, np.nan, np.nan ], - [ 1.38, 1.81, 1.38 ], - [ 1.30, 1.55, 1.15 ], + [ 8.0, 12.5, 9.0, np.nan, 10.0, 9.5 ], + [ 0.32, 0.46, 0.34, np.nan, 0.39, 0.38 ], + [ np.nan, np.nan, np.nan, np.nan, np.nan, np.nan ], + [ 1.38, 1.81, 1.38, np.nan, 1.56, 1.56 ], + [ 1.30, 1.55, 1.15, np.nan, 1.45, 1.35 ], ] ) self.assertTrue( np.array_equal( np.round( obtained, 2 ), np.round( expected, 2 ), equal_nan=True ) ) diff --git a/geos-posp/pyproject.toml b/geos-posp/pyproject.toml index 2b07d0de..ddd3d062 100644 --- a/geos-posp/pyproject.toml +++ b/geos-posp/pyproject.toml @@ -39,6 +39,7 @@ dependencies = [ "geos-geomechanics", "geos-utils", "geos-mesh", + "geos-processing", "vtk >= 9.3", "numpy >= 2.2", "pandas >= 2.2", diff --git a/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py b/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py deleted file mode 100644 index d1e44c26..00000000 --- a/geos-posp/src/PVplugins/PVGeomechanicsAnalysis.py +++ /dev/null @@ -1,365 +0,0 @@ -# SPDX-License-Identifier: Apache-2.0 -# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Martin Lemay -# ruff: noqa: E402 # disable Module level import not at top of file -import os -import sys -from typing import Union - -import numpy as np -from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] - VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, -) -from typing_extensions import Self -from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector -from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid - -dir_path = os.path.dirname( os.path.realpath( __file__ ) ) -parent_dir_path = os.path.dirname( dir_path ) -if parent_dir_path not in sys.path: - sys.path.append( parent_dir_path ) - -import PVplugins # noqa: F401 - -from geos.utils.Logger import Logger, getLogger -from geos.utils.PhysicalConstants import ( - DEFAULT_FRICTION_ANGLE_DEG, - DEFAULT_GRAIN_BULK_MODULUS, - DEFAULT_ROCK_COHESION, - WATER_DENSITY, -) -from geos_posp.filters.GeomechanicsCalculator import GeomechanicsCalculator - -__doc__ = """ -PVGeomechanicsAnalysis is a Paraview plugin that allows to compute -additional geomechanical attributes from the input mesh. - -Input and output types are vtkMultiBlockDataSet. - -To use it: - -* Load the module in Paraview: Tools>Manage Plugins...>Load new>PVGeomechanicsAnalysis. -* Select any pipeline child of the first ouput from PVExtractMergeBlocksVolume* filter. -* Search and Apply PVGeomechanicsAnalysis Filter. - -""" - - -@smproxy.filter( name="PVGeomechanicsAnalysis", label="Geos Geomechanics Analysis" ) -@smhint.xml( """""" ) -@smproperty.input( name="Input", port_index=0 ) -@smdomain.datatype( dataTypes=[ "vtkUnstructuredGrid", "vtkPointSet" ], composite_data_supported=True ) -class PVGeomechanicsAnalysis( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """Paraview plugin to compute additional geomechanical outputs. - - Input is either a vtkUnstructuredGrid or vtkPointSet with Geos - geomechanical properties. - """ - super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkPointSet" ) - - # outputs and additional parameters - self.m_computeAdvancedOutputs: bool = False - self.m_grainBulkModulus: float = DEFAULT_GRAIN_BULK_MODULUS - self.m_specificDensity: float = WATER_DENSITY - self.m_rockCohesion: float = DEFAULT_ROCK_COHESION - self.m_frictionAngle: float = DEFAULT_FRICTION_ANGLE_DEG - - # set m_logger - self.m_logger: Logger = getLogger( "Geomechanics Analysis Filter" ) - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set filter logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - - @smproperty.doublevector( - name="GrainBulkModulus", - label="Grain bulk modulus (Pa)", - default_values=DEFAULT_GRAIN_BULK_MODULUS, - panel_visibility="default", - ) - @smdomain.xml( """ - - Reference grain bulk modulus to compute Biot coefficient. - The unit is Pa. Default is Quartz bulk modulus (i.e., 38GPa). - - """ ) - def b01SetGrainBulkModulus( self: Self, value: float ) -> None: - """Set grain bulk modulus. - - Args: - value (float): grain bulk modulus (Pa) - """ - self.m_grainBulkModulus = value - self.Modified() - - def getGrainBulkModulus( self: Self ) -> float: - """Access to the grain bulk modulus value. - - Returns: - float: self.m_grainBulkModulus. - """ - return self.m_grainBulkModulus - - @smproperty.doublevector( - name="SpecificDensity", - label="Specific Density (kg/m3)", - default_values=WATER_DENSITY, - panel_visibility="default", - ) - @smdomain.xml( """ - - Reference density to compute specific gravity. - The unit is kg/m3. Default is fresh water density (i.e., 1000 kg/m3). - - """ ) - def b02SetSpecificDensity( self: Self, value: float ) -> None: - """Set specific density. - - Args: - value (float): Reference specific density (kg/m3) - """ - self.m_specificDensity = value - self.Modified() - - def getSpecificDensity( self: Self ) -> float: - """Access the specific density value. - - Returns: - float: self.m_specificDensity. - """ - return self.m_specificDensity - - @smproperty.xml( """ - - - - - """ ) - def b09GroupBasicOutputParameters( self: Self ) -> None: - """Organize groups.""" - self.Modified() - - @smproperty.intvector( - name="AdvancedOutputsUse", - label="Compute advanced geomechanical outputs", - default_values=0, - panel_visibility="default", - ) - @smdomain.xml( """ - - - Check to compute advanced geomechanical outputs including - reservoir stress paths and fracture indexes. - - """ ) - def c01SetAdvancedOutputs( self: Self, boolean: bool ) -> None: - """Set advanced output calculation option. - - Args: - boolean (bool): if True advanced outputs are computed. - """ - self.m_computeAdvancedOutputs = boolean - self.Modified() - - def getComputeAdvancedOutputs( self: Self ) -> float: - """Access the advanced outputs option. - - Returns: - float: self.m_computeAdvancedOutputs. - """ - return self.m_computeAdvancedOutputs - - @smproperty.xml( """ - - panel_visibility="default"> - - """ ) - def c09GroupAdvancedOutputsUse( self: Self ) -> None: - """Organize groups.""" - self.Modified() - - @smproperty.doublevector( - name="RockCohesion", - label="Rock Cohesion (Pa)", - default_values=DEFAULT_ROCK_COHESION, - panel_visibility="default", - ) - @smdomain.xml( """ - - Reference rock cohesion to compute critical pore pressure. - The unit is Pa.Default is fractured case (i.e., 0. Pa). - - """ ) - def d01SetRockCohesion( self: Self, value: float ) -> None: - """Set rock cohesion. - - Args: - value (float): rock cohesion (Pa) - """ - self.m_rockCohesion = value - self.Modified() - - def getRockCohesion( self: Self ) -> float: - """Get rock cohesion. - - Returns: - float: rock cohesion. - """ - return self.m_rockCohesion - - @smproperty.doublevector( - name="FrictionAngle", - label="Friction Angle (°)", - default_values=DEFAULT_FRICTION_ANGLE_DEG, - panel_visibility="default", - ) - @smdomain.xml( """ - - Reference friction angle to compute critical pore pressure. - The unit is °. Default is no friction case (i.e., 0.°). - - """ ) - def d02SetFrictionAngle( self: Self, value: float ) -> None: - """Set frition angle. - - Args: - value (float): friction angle (°) - """ - self.m_frictionAngle = value - self.Modified() - - def getFrictionAngle( self: Self ) -> float: - """Get friction angle in radian. - - Returns: - float: friction angle. - """ - return self.m_frictionAngle * np.pi / 180.0 - - @smproperty.xml( """ - - - - - - - - """ ) - def d09GroupAdvancedOutputParameters( self: Self ) -> None: - """Organize groups.""" - self.Modified() - - def RequestInformation( - self: Self, - request: vtkInformation, # noqa: F841 - inInfoVec: list[ vtkInformationVector ], # noqa: F841 - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestInformation. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - executive = self.GetExecutive() # noqa: F841 - outInfo = outInfoVec.GetInformationObject( 0 ) # noqa: F841 - return 1 - - def RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData = self.GetInputData( inInfoVec, 0, 0 ) - outData = self.GetOutputData( outInfoVec, 0 ) - assert inData is not None - if outData is None or ( not outData.IsA( inData.GetClassName() ) ): - outData = inData.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] - - def RequestData( - self: Self, - request: vtkInformation, # noqa: F841 - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestData. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - try: - self.m_logger.info( f"Apply filter {__name__}" ) - - inData = self.GetInputData( inInfoVec, 0, 0 ) - assert inData is not None - - input: Union[ vtkPointSet, vtkUnstructuredGrid ] - output: Union[ vtkPointSet, vtkUnstructuredGrid ] - if inData.IsA( "vtkUnstructuredGrid" ): - input = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) - output = vtkUnstructuredGrid.GetData( outInfoVec ) - elif inData.IsA( "vtkPointSet" ): - input = vtkPointSet.GetData( inInfoVec[ 0 ] ) - output = vtkPointSet.GetData( outInfoVec ) - else: - raise TypeError( "Error type" ) - - assert input is not None, "Input object is null" - assert output is not None, "Output object is null" - - # create new properties - geomechanicsCalculatorFilter: GeomechanicsCalculator = ( GeomechanicsCalculator() ) - geomechanicsCalculatorFilter.SetLogger( self.m_logger ) - geomechanicsCalculatorFilter.SetInputDataObject( input ) - if self.m_computeAdvancedOutputs: - geomechanicsCalculatorFilter.computeAdvancedOutputsOn() - else: - geomechanicsCalculatorFilter.computeAdvancedOutputsOff() - geomechanicsCalculatorFilter.SetGrainBulkModulus( self.getGrainBulkModulus() ) - geomechanicsCalculatorFilter.SetSpecificDensity( self.getSpecificDensity() ) - geomechanicsCalculatorFilter.SetRockCohesion( self.getRockCohesion() ) - geomechanicsCalculatorFilter.SetFrictionAngle( self.getFrictionAngle() ) - geomechanicsCalculatorFilter.Update() - output.ShallowCopy( geomechanicsCalculatorFilter.GetOutputDataObject( 0 ) ) - output.Modified() - - except AssertionError as e: - self.m_logger.error( f"{__name__} filter execution failed due to:" ) - self.m_logger.error( e ) - except Exception as e: - self.m_logger.critical( f"{__name__} filter execution failed due to:" ) - self.m_logger.critical( e, exc_info=True ) - return 1 diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py index c50e5c22..b077fe7b 100644 --- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py +++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolume.py @@ -29,9 +29,12 @@ from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, ) +from paraview.detail.loghandler import ( # type: ignore[import-not-found] + VTKHandler, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from PVplugins.PVExtractMergeBlocksVolume import PVExtractMergeBlocksVolume -from PVplugins.PVGeomechanicsAnalysis import PVGeomechanicsAnalysis +from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator __doc__ = """ PVGeomechanicsWorkflowVolume is a Paraview plugin that execute multiple filters @@ -369,15 +372,16 @@ def computeAdditionalOutputsVolume( self: Self ) -> bool: Returns: bool: True if calculation successfully eneded, False otherwise. """ - filter = PVGeomechanicsAnalysis() - filter.SetInputDataObject( self.m_volumeMesh ) - filter.b01SetGrainBulkModulus( self.getGrainBulkModulus() ) - filter.b02SetSpecificDensity( self.getSpecificDensity() ) - filter.d01SetRockCohesion( self.getRockCohesion() ) - filter.d02SetFrictionAngle( self.getFrictionAngle() ) - filter.c01SetAdvancedOutputs( self.m_computeAdvancedOutputs ) - filter.SetLogger( self.m_logger ) - filter.Update() - self.m_volumeMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + filter = GeomechanicsCalculator( self.m_volumeMesh, + computeAdvancedOutputs=self.getComputeAdvancedOutputs(), + speHandler=True ) + if not filter.logger.hasHandlers(): + filter.setLoggerHandler( VTKHandler() ) + filter.physicalConstants.grainBulkModulus = self.grainBulkModulus + filter.physicalConstants.specificDensity = self.specificDensity + filter.physicalConstants.rockCohesion = self.rockCohesion + filter.physicalConstants.frictionAngle = self.frictionAngle + filter.applyFilter() + self.m_volumeMesh.ShallowCopy( filter.getOutput() ) self.m_volumeMesh.Modified() return True diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py index 0bfae2ee..8b559688 100644 --- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py +++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurface.py @@ -21,6 +21,9 @@ from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, ) +from paraview.detail.loghandler import ( # type: ignore[import-not-found] + VTKHandler, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from geos.utils.Logger import Logger, getLogger from geos.utils.PhysicalConstants import ( @@ -32,7 +35,7 @@ ) from PVplugins.PVExtractMergeBlocksVolumeSurface import ( PVExtractMergeBlocksVolumeSurface, ) -from PVplugins.PVGeomechanicsAnalysis import PVGeomechanicsAnalysis +from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator from PVplugins.PVSurfaceGeomechanics import PVSurfaceGeomechanics __doc__ = """ @@ -370,16 +373,17 @@ def computeAdditionalOutputsVolume( self: Self ) -> bool: Returns: bool: True if calculation successfully eneded, False otherwise. """ - filter = PVGeomechanicsAnalysis() - filter.SetInputDataObject( self.m_volumeMesh ) - filter.b01SetGrainBulkModulus( self.getGrainBulkModulus() ) - filter.b02SetSpecificDensity( self.getSpecificDensity() ) - filter.d01SetRockCohesion( self.getRockCohesion() ) - filter.d02SetFrictionAngle( self.getFrictionAngle() ) - filter.c01SetAdvancedOutputs( self.m_computeAdvancedOutputs ) - filter.SetLogger( self.m_logger ) - filter.Update() - self.m_volumeMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + filter = GeomechanicsCalculator( self.m_volumeMesh, + computeAdvancedOutputs=self.getComputeAdvancedOutputs(), + speHandler=True ) + if not filter.logger.hasHandlers(): + filter.setLoggerHandler( VTKHandler() ) + filter.physicalConstants.grainBulkModulus = self.grainBulkModulus + filter.physicalConstants.specificDensity = self.specificDensity + filter.physicalConstants.rockCohesion = self.rockCohesion + filter.physicalConstants.frictionAngle = self.frictionAngle + filter.applyFilter() + self.m_volumeMesh.ShallowCopy( filter.getOutput() ) self.m_volumeMesh.Modified() return True diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py index fff90856..715de74c 100644 --- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py +++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeSurfaceWell.py @@ -21,6 +21,9 @@ from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, ) +from paraview.detail.loghandler import ( # type: ignore[import-not-found] + VTKHandler, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from geos.utils.Logger import Logger, getLogger from geos.utils.PhysicalConstants import ( @@ -32,7 +35,7 @@ ) from PVplugins.PVExtractMergeBlocksVolumeSurfaceWell import ( PVExtractMergeBlocksVolumeSurfaceWell, ) -from PVplugins.PVGeomechanicsAnalysis import PVGeomechanicsAnalysis +from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator from PVplugins.PVSurfaceGeomechanics import PVSurfaceGeomechanics __doc__ = """ @@ -379,16 +382,17 @@ def computeAdditionalOutputsVolume( self: Self ) -> bool: Returns: bool: True if calculation successfully eneded, False otherwise. """ - filter = PVGeomechanicsAnalysis() - filter.SetInputDataObject( self.m_volumeMesh ) - filter.b01SetGrainBulkModulus( self.getGrainBulkModulus() ) - filter.b02SetSpecificDensity( self.getSpecificDensity() ) - filter.d01SetRockCohesion( self.getRockCohesion() ) - filter.d02SetFrictionAngle( self.getFrictionAngle() ) - filter.c01SetAdvancedOutputs( self.m_computeAdvancedOutputs ) - filter.SetLogger( self.m_logger ) - filter.Update() - self.m_volumeMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + filter = GeomechanicsCalculator( self.m_volumeMesh, + computeAdvancedOutputs=self.getComputeAdvancedOutputs(), + speHandler=True ) + if not filter.logger.hasHandlers(): + filter.setLoggerHandler( VTKHandler() ) + filter.physicalConstants.grainBulkModulus = self.grainBulkModulus + filter.physicalConstants.specificDensity = self.specificDensity + filter.physicalConstants.rockCohesion = self.rockCohesion + filter.physicalConstants.frictionAngle = self.frictionAngle + filter.applyFilter() + self.m_volumeMesh.ShallowCopy( filter.getOutput() ) self.m_volumeMesh.Modified() return True diff --git a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py index 210fd933..0d5829dd 100644 --- a/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py +++ b/geos-posp/src/PVplugins/PVGeomechanicsWorkflowVolumeWell.py @@ -21,6 +21,9 @@ from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, ) +from paraview.detail.loghandler import ( # type: ignore[import-not-found] + VTKHandler, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py from geos.utils.Logger import Logger, getLogger from geos.utils.PhysicalConstants import ( @@ -33,7 +36,7 @@ from PVplugins.PVExtractMergeBlocksVolumeWell import ( PVExtractMergeBlocksVolumeWell, ) -from PVplugins.PVGeomechanicsAnalysis import PVGeomechanicsAnalysis +from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator __doc__ = """ PVGeomechanicsWorkflowVolumeWell is a Paraview plugin that execute @@ -384,15 +387,16 @@ def computeAdditionalOutputsVolume( self: Self ) -> bool: Returns: bool: True if calculation successfully eneded, False otherwise. """ - filter = PVGeomechanicsAnalysis() - filter.SetInputDataObject( self.m_volumeMesh ) - filter.b01SetGrainBulkModulus( self.getGrainBulkModulus() ) - filter.b02SetSpecificDensity( self.getSpecificDensity() ) - filter.d01SetRockCohesion( self.getRockCohesion() ) - filter.d02SetFrictionAngle( self.getFrictionAngle() ) - filter.c01SetAdvancedOutputs( self.m_computeAdvancedOutputs ) - filter.SetLogger( self.m_logger ) - filter.Update() - self.m_volumeMesh.ShallowCopy( filter.GetOutputDataObject( 0 ) ) + filter = GeomechanicsCalculator( self.m_volumeMesh, + computeAdvancedOutputs=self.getComputeAdvancedOutputs(), + speHandler=True ) + if not filter.logger.hasHandlers(): + filter.setLoggerHandler( VTKHandler() ) + filter.physicalConstants.grainBulkModulus = self.grainBulkModulus + filter.physicalConstants.specificDensity = self.specificDensity + filter.physicalConstants.rockCohesion = self.rockCohesion + filter.physicalConstants.frictionAngle = self.frictionAngle + filter.applyFilter() + self.m_volumeMesh.ShallowCopy( filter.getOutput() ) self.m_volumeMesh.Modified() return True diff --git a/geos-posp/src/PVplugins/__init__.py b/geos-posp/src/PVplugins/__init__.py index 6271e4ad..37e7c3cf 100644 --- a/geos-posp/src/PVplugins/__init__.py +++ b/geos-posp/src/PVplugins/__init__.py @@ -5,7 +5,7 @@ dir_path = os.path.dirname( os.path.realpath( __file__ ) ) python_root = '../../..' -python_modules = ( 'geos-posp', 'geos-utils', 'geos-geomechanics', 'geos-mesh' ) +python_modules = ( 'geos-posp', 'geos-utils', 'geos-geomechanics', 'geos-mesh', 'geos-processing' ) for m in python_modules: m_path = os.path.abspath( os.path.join( dir_path, python_root, m, 'src' ) ) diff --git a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py b/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py deleted file mode 100644 index fb80122a..00000000 --- a/geos-posp/src/geos_posp/filters/GeomechanicsCalculator.py +++ /dev/null @@ -1,1382 +0,0 @@ -# SPDX-License-Identifier: Apache-2.0 -# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Martin Lemay -# ruff: noqa: E402 # disable Module level import not at top of file -from typing import Union - -import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts -import numpy as np -import numpy.typing as npt -from geos.utils.GeosOutputsConstants import ( - AttributeEnum, - ComponentNameEnum, - GeosMeshOutputsEnum, - PostProcessingOutputsEnum, -) -from geos.utils.Logger import Logger, getLogger -from geos.utils.PhysicalConstants import ( - DEFAULT_FRICTION_ANGLE_RAD, - DEFAULT_GRAIN_BULK_MODULUS, - DEFAULT_ROCK_COHESION, - GRAVITY, - WATER_DENSITY, -) -from typing_extensions import Self -from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase -from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector -from vtkmodules.vtkCommonDataModel import ( - vtkDataSet, - vtkPointSet, - vtkUnstructuredGrid, -) -from vtkmodules.vtkFiltersCore import vtkCellCenters -from geos.mesh.utils.arrayModifiers import createAttribute -from geos.mesh.utils.arrayHelpers import ( - getArrayInObject, - getComponentNames, - isAttributeInObject, -) - -__doc__ = """ -GeomechanicsCalculator module is a vtk filter that allows to compute additional -Geomechanical properties from existing ones. - -GeomechanicsCalculator filter inputs are either vtkPointSet or vtkUnstructuredGrid -and returned object is of same type as input. - -To use the filter: - -.. code-block:: python - - import numpy as np - from geos_posp.filters.GeomechanicsCalculator import GeomechanicsCalculator - - # filter inputs - logger :Logger - # input object - input :Union[vtkPointSet, vtkUnstructuredGrid] - # grain bulk modulus in Pa (e.g., use a very high value to get a Biot coefficient equal to 1) - grainBulkModulus :float = 1e26 - # Reference density to compute specific gravity (e.g. fresh water density) in kg.m^-3 - specificDensity :float = 1000. - # rock cohesion in Pa - rockCohesion :float = 1e8 - # friction angle in ° - frictionAngle :float = 10 * np.pi / 180. - - # instanciate the filter - geomechanicsCalculatorFilter :GeomechanicsCalculator = GeomechanicsCalculator() - - # set filter attributes - # set logger - geomechanicsCalculatorFilter.SetLogger(logger) - # set input object - geomechanicsCalculatorFilter.SetInputDataObject(input) - # set computeAdvancedOutputsOn or computeAdvancedOutputsOff to compute or - # not advanced outputs... - geomechanicsCalculatorFilter.computeAdvancedOutputsOn() - # set oter parameters - geomechanicsCalculatorFilter.SetGrainBulkModulus(grainBulkModulus) - geomechanicsCalculatorFilter.SetSpecificDensity(specificDensity) - # rock cohesion and friction angle are used for advanced outputs only - geomechanicsCalculatorFilter.SetRockCohesion(rockCohesion) - geomechanicsCalculatorFilter.SetFrictionAngle(frictionAngle) - # do calculations - geomechanicsCalculatorFilter.Update() - # get filter output (same type as input) - output :Union[vtkPointSet, vtkUnstructuredGrid] = geomechanicsCalculatorFilter.GetOutputDataObject(0) -""" - -TYPE_ERROR_MESSAGE = ( "Input object must by either a vtkPointSet or a vtkUntructuredGrid." ) -UNDEFINED_ATTRIBUTE_MESSAGE = " attribute is undefined." - - -class GeomechanicsCalculator( VTKPythonAlgorithmBase ): - - def __init__( self: Self ) -> None: - """VTK Filter to perform Geomechanical output computation. - - Input object is either a vtkPointSet or a vtkUntructuredGrid. - """ - super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkDataSet" ) # type: ignore[no-untyped-call] - - self.m_output: Union[ vtkPointSet, vtkUnstructuredGrid ] - - # additional parameters - self.m_computeAdvancedOutputs: bool = False - self.m_grainBulkModulus: float = DEFAULT_GRAIN_BULK_MODULUS - self.m_specificDensity: float = WATER_DENSITY - self.m_rockCohesion: float = DEFAULT_ROCK_COHESION - self.m_frictionAngle: float = DEFAULT_FRICTION_ANGLE_RAD - - # computation results - self.m_elasticModuliComputed: bool = False - self.m_biotCoefficientComputed: bool = False - self.m_compressibilityComputed: bool = False - self.m_effectiveStressComputed: bool = False - self.m_totalStressComputed: bool = False - self.m_effectiveStressRatioOedComputed: bool = False - - # will compute resuls if m_ready is True (updated by initFilter method) - self.m_ready: bool = False - # attributes are either on points or on cells - self.m_attributeOnPoints: bool = False - # elastic moduli are either bulk and Shear moduli (m_computeYoungPoisson=True) - # or young Modulus and poisson's ratio (m_computeYoungPoisson=False) - self.m_computeYoungPoisson: bool = True - - # set m_logger - self.m_logger: Logger = getLogger( "Geomechanics calculator" ) - - def FillInputPortInformation( self: Self, port: int, info: vtkInformation ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestInformation. - - Args: - port (int): input port - info (vtkInformationVector): info - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - if port == 0: - info.Set( self.INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet" ) - return 1 - - def RequestInformation( - self: Self, - request: vtkInformation, # noqa: F841 - inInfoVec: list[ vtkInformationVector ], # noqa: F841 - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestInformation. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - executive = self.GetExecutive() # noqa: F841 - outInfo = outInfoVec.GetInformationObject( 0 ) # noqa: F841 - return 1 - - def RequestDataObject( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestDataObject. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - inData = self.GetInputData( inInfoVec, 0, 0 ) # type: ignore[no-untyped-call] - outData = self.GetOutputData( outInfoVec, 0 ) # type: ignore[no-untyped-call] - assert inData is not None - if outData is None or ( not outData.IsA( inData.GetClassName() ) ): - outData = inData.NewInstance() - outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) - return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore - - def RequestData( - self: Self, - request: vtkInformation, - inInfoVec: list[ vtkInformationVector ], - outInfoVec: vtkInformationVector, - ) -> int: - """Inherited from VTKPythonAlgorithmBase::RequestData. - - Args: - request (vtkInformation): request - inInfoVec (list[vtkInformationVector]): input objects - outInfoVec (vtkInformationVector): output objects - - Returns: - int: 1 if calculation successfully ended, 0 otherwise. - """ - try: - input: vtkDataSet = vtkDataSet.GetData( inInfoVec[ 0 ] ) - assert input is not None, "Input object is null." - - # initialize output objects - self.m_output = self.GetOutputData( outInfoVec, 0 ) # type: ignore[no-untyped-call] - assert self.m_output is not None, "Output object is null." - self.m_output.ShallowCopy( input ) - - # check the input and update self.m_ready, m_attributeOnPoints and m_computeBulkAndShear - self.initFilter() - if not self.m_ready: - mess: str = ( "Mandatory properties are missing to compute geomechanical outputs:" ) - mess += ( f"Either {PostProcessingOutputsEnum.YOUNG_MODULUS.attributeName} " - f"and {PostProcessingOutputsEnum.POISSON_RATIO.attributeName} or " - f"{GeosMeshOutputsEnum.BULK_MODULUS.attributeName} and " - f"{GeosMeshOutputsEnum.SHEAR_MODULUS.attributeName} must be " - f"present in the data as well as the " - f"{GeosMeshOutputsEnum.STRESS_EFFECTIVE.attributeName} attribute." ) - self.m_logger.error( mess ) - return 0 - - self.computeBasicOutputs() - if self.m_computeAdvancedOutputs: - self.computeAdvancedOutputs() - - except AssertionError as e: - mess1: str = "Geomechanical attribute calculation failed due to:" - self.m_logger.error( mess1 ) - self.m_logger.error( str( e ) ) - return 0 - except Exception as e: - mess2: str = "Geomechanical attribut calculation failed due to:" - self.m_logger.critical( mess2 ) - self.m_logger.critical( e, exc_info=True ) - return 0 - return 1 - - def SetLogger( self: Self, logger: Logger ) -> None: - """Set the logger. - - Args: - logger (Logger): logger - """ - self.m_logger = logger - self.Modified() - - def computeAdvancedOutputsOn( self: Self ) -> None: - """Activate advanced outputs calculation.""" - self.m_computeAdvancedOutputs = True - self.Modified() - - def computeAdvancedOutputsOff( self: Self ) -> None: - """Deactivate advanced outputs calculation.""" - self.m_computeAdvancedOutputs = False - self.Modified() - - def SetGrainBulkModulus( self: Self, grainBulkModulus: float ) -> None: - """Set the grain bulk modulus. - - Args: - grainBulkModulus (float): grain bulk modulus - """ - self.m_grainBulkModulus = grainBulkModulus - self.Modified() - - def SetSpecificDensity( self: Self, specificDensity: float ) -> None: - """Set the specific density. - - Args: - specificDensity (float): pecific density - """ - self.m_specificDensity = specificDensity - self.Modified() - - def SetRockCohesion( self: Self, rockCohesion: float ) -> None: - """Set the rock cohesion. - - Args: - rockCohesion (float): rock cohesion - """ - self.m_rockCohesion = rockCohesion - self.Modified() - - def SetFrictionAngle( self: Self, frictionAngle: float ) -> None: - """Set the friction angle. - - Args: - frictionAngle (float): friction angle (rad) - """ - self.m_frictionAngle = frictionAngle - self.Modified() - - def getOutputType( self: Self ) -> str: - """Get output object type. - - Returns: - str: type of output object. - """ - output: vtkDataSet = self.GetOutputDataObject( 0 ) - assert output is not None, "Output is null." - return output.GetClassName() - - def resetDefaultValues( self: Self ) -> None: - """Reset filter parameters to the default values.""" - self.m_computeAdvancedOutputs = False - self.m_grainBulkModulus = DEFAULT_GRAIN_BULK_MODULUS - self.m_specificDensity = WATER_DENSITY - self.m_rockCohesion = DEFAULT_ROCK_COHESION - self.m_frictionAngle = DEFAULT_FRICTION_ANGLE_RAD - - self.m_elasticModuliComputed = False - self.m_biotCoefficientComputed = False - self.m_compressibilityComputed = False - self.m_effectiveStressComputed = False - self.m_totalStressComputed = False - self.m_effectiveStressRatioOedComputed = False - self.m_ready = False - self.Modified() - - def initFilter( self: Self ) -> None: - """Check that mandatory attributes are present in the data set. - - Determine if attributes are on cells or on Points. - Set self.m_ready = True if all data is ok, False otherwise - """ - # check attributes are on cells, or on points otherwise - attributeOnPoints: bool = False - attributeOnCells: bool = self.checkMandatoryAttributes( False ) - if not attributeOnCells: - attributeOnPoints = self.checkMandatoryAttributes( True ) - - self.m_ready = attributeOnPoints or attributeOnCells - self.m_attributeOnPoints = attributeOnPoints - - def checkMandatoryAttributes( self: Self, onPoints: bool ) -> bool: - """Check that mandatory attributes are present in the mesh. - - The mesh must contains either the young Modulus and Poisson's ratio - (m_computeYoungPoisson=False) or the bulk and shear moduli - (m_computeYoungPoisson=True) - - Args: - onPoints (bool): attributes are on points (True) or on cells (False) - - Returns: - bool: True if all needed attributes are present, False otherwise - """ - youngModulusAttributeName: str = ( PostProcessingOutputsEnum.YOUNG_MODULUS.attributeName ) - poissonRatioAttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO.attributeName ) - bulkModulusAttributeName: str = GeosMeshOutputsEnum.BULK_MODULUS.attributeName - shearModulusAttributeName: str = GeosMeshOutputsEnum.SHEAR_MODULUS.attributeName - effectiveStressAttributeName: str = ( GeosMeshOutputsEnum.STRESS_EFFECTIVE.attributeName ) - - self.m_computeYoungPoisson = not isAttributeInObject( self.m_output, youngModulusAttributeName, - onPoints ) or not isAttributeInObject( - self.m_output, poissonRatioAttributeName, onPoints ) - - # if none of elastic moduli is present, return False - if self.m_computeYoungPoisson and ( - not isAttributeInObject( self.m_output, bulkModulusAttributeName, onPoints ) - or not isAttributeInObject( self.m_output, shearModulusAttributeName, onPoints ) ): - return False - - # check effective Stress is present - ret: bool = isAttributeInObject( self.m_output, effectiveStressAttributeName, onPoints ) - return ret - - def computeBasicOutputs( self: Self ) -> bool: - """Compute basic geomechanical outputs. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - if not self.m_ready: - return False - - try: - self.m_elasticModuliComputed = self.computeElasticModulus() - if not self.m_elasticModuliComputed: - mess: str = ( "Geomechanical outputs cannot be computed without elastic moduli." ) - self.m_logger.error( mess ) - return False - - self.m_biotCoefficientComputed = self.computeBiotCoefficient() - if not self.m_biotCoefficientComputed: - mess2: str = ( "Total stress, elastic strain, and advanced geomechanical " + - "outputs cannot be computed without Biot coefficient." ) - self.m_logger.warning( mess2 ) - - self.m_compressibilityComputed = self.computeCompressibilityCoefficient() - - self.m_effectiveStressComputed = self.computeRealEffectiveStressRatio() - if not self.m_effectiveStressComputed: - mess3: str = ( "Total stress, elastic strain, and advanced geomechanical " + - "outputs cannot be computed without effective stress." ) - self.m_logger.warning( mess3 ) - - specificGravityComputed: bool = self.computeSpecificGravity() - - # TODO: deactivate lithostatic stress calculation until right formula - litostaticStressComputed: bool = True # self.computeLitostaticStress() - - elasticStrainComputed: bool = False - if self.m_effectiveStressComputed: - if self.m_biotCoefficientComputed: - self.m_totalStressComputed = self.computeTotalStresses() - if self.m_elasticModuliComputed: - elasticStrainComputed = self.computeElasticStrain() - - reservoirStressPathOedComputed: bool = False - if self.m_elasticModuliComputed: - # oedometric DRSP (effective stress ratio in oedometric conditions) - self.m_effectiveStressRatioOedComputed = ( self.computeEffectiveStressRatioOed() ) - - if self.m_biotCoefficientComputed: - reservoirStressPathOedComputed = ( self.computeReservoirStressPathOed() ) - - reservoirStressPathRealComputed: bool = False - if self.m_totalStressComputed: - reservoirStressPathRealComputed = self.computeReservoirStressPathReal() - - if ( self.m_elasticModuliComputed and self.m_biotCoefficientComputed and self.m_compressibilityComputed - and self.m_effectiveStressComputed and specificGravityComputed and elasticStrainComputed - and litostaticStressComputed and self.m_totalStressComputed and self.m_effectiveStressRatioOedComputed - and reservoirStressPathRealComputed and reservoirStressPathRealComputed - and reservoirStressPathOedComputed and reservoirStressPathRealComputed ): - mess4: str = ( "All geomechanical basic outputs were successfully computed." ) - self.m_logger.info( mess4 ) - else: - mess5: str = "Some geomechanical basic outputs were not computed." - self.m_logger.warning( mess5 ) - - except AssertionError as e: - mess6: str = ( "Some of the geomechanical basic outputs were " + "not computed due to:" ) - self.m_logger.error( mess6 ) - self.m_logger.error( str( e ) ) - return False - except Exception as e: - mess7: str = ( "Some of the geomechanical basic outputs were " + "not computed due to:" ) - self.m_logger.critical( mess7 ) - self.m_logger.critical( e, exc_info=True ) - return False - return True - - def computeAdvancedOutputs( self: Self ) -> bool: - """Compute advanced geomechanical outputs. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - if not self.m_ready: - return False - - try: - fractureIndexesComputed: bool = False - criticalPorePressure: bool = False - if self.m_totalStressComputed: - fractureIndexesComputed = self.computeCriticalTotalStressRatio() - criticalPorePressure = self.computeCriticalPorePressure() - - if ( self.m_effectiveStressRatioOedComputed and fractureIndexesComputed and criticalPorePressure ): - mess: str = ( "All geomechanical advanced outputs were " + "successfully computed." ) - self.m_logger.info( mess ) - else: - mess0: str = ( "Some geomechanical advanced outputs were " + "not computed." ) - self.m_logger.warning( mess0 ) - - except AssertionError as e: - mess1: str = ( "Some of the geomechanical basic outputs were " + "not computed due to:" ) - self.m_logger.error( mess1 ) - self.m_logger.error( str( e ) ) - return False - except Exception as e: - mess2: str = ( "Some of the geomechanical advanced outputs " + "were not computed due to:" ) - self.m_logger.critical( mess2 ) - self.m_logger.critical( e, exc_info=True ) - return False - return True - - def computeElasticModulus( self: Self ) -> bool: - """Compute elastic moduli. - - Young modulus and the poisson ratio computed from shear and bulk moduli - if needed. - - Returns: - bool: True if elastic moduli are already present or if calculation - successfully ended, False otherwise. - """ - if self.m_computeYoungPoisson: - return self.computeElasticModulusFromBulkShear() - return self.computeElasticModulusFromYoungPoisson() - - def computeElasticModulusFromBulkShear( self: Self ) -> bool: - """Compute Young modulus Poisson's ratio from bulk and shear moduli. - - Returns: - bool: True if calculation successfully ended, False otherwise - """ - youngModulusAttributeName: str = ( PostProcessingOutputsEnum.YOUNG_MODULUS.attributeName ) - poissonRatioAttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO.attributeName ) - bulkModulusAttributeName: str = GeosMeshOutputsEnum.BULK_MODULUS.attributeName - shearModulusAttributeName: str = GeosMeshOutputsEnum.SHEAR_MODULUS.attributeName - - ret: bool = True - bulkModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, bulkModulusAttributeName, - self.m_attributeOnPoints ) - - shearModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, shearModulusAttributeName, - self.m_attributeOnPoints ) - try: - assert bulkModulus is not None, ( f"{bulkModulusAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert shearModulus is not None, ( f"{shearModulusAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - except AssertionError as e: - self.m_logger.error( "Elastic moduli were not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - try: - youngModulus: npt.NDArray[ np.float64 ] = fcts.youngModulus( bulkModulus, shearModulus ) - # assert np.any(youngModulus < 0), ("Young modulus yields negative " + - # "values. Check Bulk and Shear modulus values.") - createAttribute( - self.m_output, - youngModulus, - youngModulusAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Young modulus was not computed due to:" ) - self.m_logger.error( str( e ) ) - ret = False - - try: - poissonRatio: npt.NDArray[ np.float64 ] = fcts.poissonRatio( bulkModulus, shearModulus ) - # assert np.any(poissonRatio < 0), ("Poisson ratio yields negative " + - # "values. Check Bulk and Shear modulus values.") - createAttribute( - self.m_output, - poissonRatio, - poissonRatioAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Poisson's ratio was not computed due to:" ) - self.m_logger.error( str( e ) ) - ret = False - - return ret - - def computeElasticModulusFromYoungPoisson( self: Self ) -> bool: - """Compute bulk modulus from Young Modulus and Poisson's ratio. - - Returns: - bool: True if bulk modulus was wuccessfully computed, False otherwise - """ - try: - youngModulusAttributeName: str = ( PostProcessingOutputsEnum.YOUNG_MODULUS.attributeName ) - poissonRatioAttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO.attributeName ) - bulkModulusAttributeName: str = ( GeosMeshOutputsEnum.BULK_MODULUS.attributeName ) - if not isAttributeInObject( self.m_output, bulkModulusAttributeName, self.m_attributeOnPoints ): - youngModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, youngModulusAttributeName, - self.m_attributeOnPoints ) - poissonRatio: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, poissonRatioAttributeName, - self.m_attributeOnPoints ) - - assert youngModulus is not None, ( f"{youngModulusAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert poissonRatio is not None, ( f"{poissonRatioAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - bulkModulus: npt.NDArray[ np.float64 ] = fcts.bulkModulus( youngModulus, poissonRatio ) - # assert np.any(bulkModulus < 0), ("Bulk modulus yields negative " + - # "values. Check Young modulus and Poisson ratio values.") - ret: bool = createAttribute( - self.m_output, - bulkModulus, - bulkModulusAttributeName, - (), - self.m_attributeOnPoints, - ) - return ret - - except AssertionError as e: - self.m_logger.error( "Bulk modulus was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - return True - - def computeBiotCoefficient( self: Self ) -> bool: - """Compute Biot coefficient from default and grain bulk modulus. - - Returns: - bool: True if calculation successfully ended, False otherwise. - """ - biotCoefficientAttributeName: str = ( PostProcessingOutputsEnum.BIOT_COEFFICIENT.attributeName ) - if not isAttributeInObject( self.m_output, biotCoefficientAttributeName, self.m_attributeOnPoints ): - bulkModulusAttributeName: str = ( GeosMeshOutputsEnum.BULK_MODULUS.attributeName ) - bulkModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, bulkModulusAttributeName, - self.m_attributeOnPoints ) - try: - assert bulkModulus is not None, ( f"{bulkModulusAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - - biotCoefficient: npt.NDArray[ np.float64 ] = fcts.biotCoefficient( self.m_grainBulkModulus, - bulkModulus ) - createAttribute( - self.m_output, - biotCoefficient, - biotCoefficientAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Biot coefficient was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeCompressibilityCoefficient( self: Self ) -> bool: - """Compute compressibility coefficient from simulation outputs. - - Compressibility coefficient is computed from Poisson's ratio, bulk - modulus, Biot coefficient and Porosity. - - Returns: - bool: True if the attribute is correctly created, False otherwise. - """ - compressibilityAttributeName: str = ( PostProcessingOutputsEnum.COMPRESSIBILITY.attributeName ) - bulkModulusAttributeName: str = GeosMeshOutputsEnum.BULK_MODULUS.attributeName - bulkModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, bulkModulusAttributeName, - self.m_attributeOnPoints ) - porosityAttributeName: str = GeosMeshOutputsEnum.POROSITY.attributeName - porosity: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, porosityAttributeName, - self.m_attributeOnPoints ) - porosityInitialAttributeName: str = ( GeosMeshOutputsEnum.POROSITY_INI.attributeName ) - porosityInitial: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, porosityInitialAttributeName, - self.m_attributeOnPoints ) - if not isAttributeInObject( self.m_output, compressibilityAttributeName, self.m_attributeOnPoints ): - poissonRatioAttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO.attributeName ) - poissonRatio: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, poissonRatioAttributeName, - self.m_attributeOnPoints ) - biotCoefficientAttributeName: str = ( PostProcessingOutputsEnum.BIOT_COEFFICIENT.attributeName ) - biotCoefficient: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, biotCoefficientAttributeName, - self.m_attributeOnPoints ) - - try: - assert poissonRatio is not None, ( f"{poissonRatioAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert bulkModulus is not None, ( f"{bulkModulusAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert biotCoefficient is not None, ( f"{biotCoefficientAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert porosity is not None, ( f"{porosityAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - - compressibility: npt.NDArray[ np.float64 ] = fcts.compressibility( poissonRatio, bulkModulus, - biotCoefficient, porosity ) - createAttribute( - self.m_output, - compressibility, - compressibilityAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Compressibility was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - # oedometric compressibility - compressibilityOedAttributeName: str = ( PostProcessingOutputsEnum.COMPRESSIBILITY_OED.attributeName ) - if not isAttributeInObject( self.m_output, compressibilityOedAttributeName, self.m_attributeOnPoints ): - shearModulusAttributeName: str = ( GeosMeshOutputsEnum.SHEAR_MODULUS.attributeName ) - shearModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, shearModulusAttributeName, - self.m_attributeOnPoints ) - try: - assert poissonRatio is not None, ( f"{poissonRatioAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert bulkModulus is not None, ( f"{bulkModulusAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert porosityInitial is not None, ( f"{porosityInitialAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - compressibilityOed: npt.NDArray[ np.float64 ] = fcts.compressibilityOed( - bulkModulus, shearModulus, porosityInitial ) - createAttribute( - self.m_output, - compressibilityOed, - compressibilityOedAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Oedometric Compressibility was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - # real compressibility - compressibilityRealAttributeName: str = ( PostProcessingOutputsEnum.COMPRESSIBILITY_REAL.attributeName ) - if not isAttributeInObject( self.m_output, compressibilityRealAttributeName, self.m_attributeOnPoints ): - deltaPressureAttributeName: str = ( GeosMeshOutputsEnum.DELTA_PRESSURE.attributeName ) - deltaPressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, deltaPressureAttributeName, - self.m_attributeOnPoints ) - try: - assert deltaPressure is not None, ( f"{deltaPressureAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert porosityInitial is not None, ( f"{porosityInitialAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert porosity is not None, ( f"{porosityAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - compressibilityReal: npt.NDArray[ np.float64 ] = fcts.compressibilityReal( - deltaPressure, porosity, porosityInitial ) - createAttribute( - self.m_output, - compressibilityReal, - compressibilityRealAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Real compressibility was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeSpecificGravity( self: Self ) -> bool: - """Create Specific gravity attribute. - - Specific gravity is computed from rock density attribute and specific - density input. - - Returns: - bool: True if the attribute is correctly created, False otherwise. - """ - densityAttributeName: str = GeosMeshOutputsEnum.ROCK_DENSITY.attributeName - density: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, densityAttributeName, - self.m_attributeOnPoints ) - specificGravityAttributeName: str = ( PostProcessingOutputsEnum.SPECIFIC_GRAVITY.attributeName ) - if not isAttributeInObject( self.m_output, specificGravityAttributeName, self.m_attributeOnPoints ): - try: - assert density is not None, ( f"{densityAttributeName} " + UNDEFINED_ATTRIBUTE_MESSAGE ) - - specificGravity: npt.NDArray[ np.float64 ] = fcts.specificGravity( density, self.m_specificDensity ) - createAttribute( - self.m_output, - specificGravity, - specificGravityAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Specific gravity was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeRealEffectiveStressRatio( self: Self ) -> bool: - """Compute effective stress ratio. - - Returns: - bool: True if calculation successfully ended, False otherwise. - """ - # test if effective stress is in the table - if not isAttributeInObject( - self.m_output, - GeosMeshOutputsEnum.STRESS_EFFECTIVE.attributeName, - self.m_attributeOnPoints, - ): - self.m_logger.error( "Effective stress is not in input data." ) - return False - - # real effective stress ratio - return self.computeStressRatioReal( - GeosMeshOutputsEnum.STRESS_EFFECTIVE, - PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_REAL, - ) - - def computeTotalStresses( self: Self ) -> bool: - """Compute total stress total stress ratio. - - Total stress is computed at the initial and current time steps. - total stress ratio is computed at current time step only. - - Returns: - bool: True if calculation successfully ended, False otherwise. - """ - # compute total stress at initial time step - totalStressT0AttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL_INITIAL.attributeName ) - if not isAttributeInObject( self.m_output, totalStressT0AttributeName, self.m_attributeOnPoints ): - self.computeTotalStressInitial() - - # compute total stress at current time step - totalStressAttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL.attributeName ) - if not isAttributeInObject( self.m_output, totalStressAttributeName, self.m_attributeOnPoints ): - try: - effectiveStressAttributeName: str = ( GeosMeshOutputsEnum.STRESS_EFFECTIVE.attributeName ) - effectiveStress: npt.NDArray[ np.float64 ] = getArrayInObject( - self.m_output, - effectiveStressAttributeName, - self.m_attributeOnPoints, - ) - assert effectiveStress is not None, ( f"{effectiveStressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - pressureAttributeName: str = GeosMeshOutputsEnum.PRESSURE.attributeName - pressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, pressureAttributeName, - self.m_attributeOnPoints ) - assert pressure is not None, ( f"{pressureAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - biotCoefficientAttributeName: str = ( PostProcessingOutputsEnum.BIOT_COEFFICIENT.attributeName ) - biotCoefficient: npt.NDArray[ np.float64 ] = getArrayInObject( - self.m_output, - biotCoefficientAttributeName, - self.m_attributeOnPoints, - ) - assert biotCoefficient is not None, ( f"{biotCoefficientAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - self.doComputeTotalStress( effectiveStress, pressure, biotCoefficient, totalStressAttributeName ) - self.computeStressRatioReal( - PostProcessingOutputsEnum.STRESS_TOTAL, - PostProcessingOutputsEnum.STRESS_TOTAL_RATIO_REAL, - ) - - except AssertionError as e: - self.m_logger.error( "Total stress at current time step was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeTotalStressInitial( self: Self ) -> bool: - """Compute total stress at initial time step. - - Returns: - bool: True if calculation successfully ended, False otherwise. - """ - try: - - # compute elastic moduli at initial time step - bulkModulusT0: npt.NDArray[ np.float64 ] - - bulkModulusT0AttributeName: str = ( PostProcessingOutputsEnum.BULK_MODULUS_INITIAL.attributeName ) - youngModulusT0AttributeName: str = ( PostProcessingOutputsEnum.YOUNG_MODULUS_INITIAL.attributeName ) - poissonRatioT0AttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO_INITIAL.attributeName ) - # get bulk modulus at initial time step - if isAttributeInObject( self.m_output, bulkModulusT0AttributeName, self.m_attributeOnPoints ): - - bulkModulusT0 = getArrayInObject( self.m_output, bulkModulusT0AttributeName, self.m_attributeOnPoints ) - # or compute it from Young and Poisson if not an attribute - elif isAttributeInObject( self.m_output, youngModulusT0AttributeName, - self.m_attributeOnPoints ) and isAttributeInObject( - self.m_output, poissonRatioT0AttributeName, self.m_attributeOnPoints ): - - youngModulusT0: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, - youngModulusT0AttributeName, - self.m_attributeOnPoints ) - poissonRatioT0: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, - poissonRatioT0AttributeName, - self.m_attributeOnPoints ) - assert poissonRatioT0 is not None, "Poisson's ratio is undefined." - assert youngModulusT0 is not None, "Young modulus is undefined." - bulkModulusT0 = fcts.bulkModulus( youngModulusT0, poissonRatioT0 ) - else: - raise AssertionError( "Elastic moduli at initial time are absent." ) - - assert ( bulkModulusT0 is not None ), "Bulk modulus at initial time step is undefined." - - # compute Biot at initial time step - biotCoefficientT0: npt.NDArray[ np.float64 ] = fcts.biotCoefficient( self.m_grainBulkModulus, - bulkModulusT0 ) - assert ( biotCoefficientT0 is not None ), "Biot coefficient at initial time step is undefined." - - # compute total stress at initial time step - # get effective stress attribute - effectiveStressAttributeName: str = ( PostProcessingOutputsEnum.STRESS_EFFECTIVE_INITIAL.attributeName ) - effectiveStress: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, effectiveStressAttributeName, - self.m_attributeOnPoints ) - assert effectiveStress is not None, ( f"{effectiveStressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - # get pressure attribute - pressureAttributeName: str = GeosMeshOutputsEnum.PRESSURE.attributeName - pressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, pressureAttributeName, - self.m_attributeOnPoints ) - pressureT0: npt.NDArray[ np.float64 ] - # case pressureT0 is None, total stress = effective stress - # (managed by doComputeTotalStress function) - if pressure is not None: - # get delta pressure to recompute pressure at initial time step (pressureTo) - deltaPressureAttributeName: str = ( GeosMeshOutputsEnum.DELTA_PRESSURE.attributeName ) - deltaPressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, deltaPressureAttributeName, - self.m_attributeOnPoints ) - assert deltaPressure is not None, ( f"{deltaPressureAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - pressureT0 = pressure - deltaPressure - - totalStressT0AttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL_INITIAL.attributeName ) - return self.doComputeTotalStress( - effectiveStress, - pressureT0, - biotCoefficientT0, - totalStressT0AttributeName, - ) - - except AssertionError as e: - self.m_logger.error( "Total stress at initial time step was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - def doComputeTotalStress( - self: Self, - effectiveStress: npt.NDArray[ np.float64 ], - pressure: Union[ npt.NDArray[ np.float64 ], None ], - biotCoefficient: npt.NDArray[ np.float64 ], - totalStressAttributeName: str, - ) -> bool: - """Compute total stress from effective stress and pressure. - - Args: - effectiveStress (npt.NDArray[np.float64]): effective stress - pressure (npt.NDArray[np.float64] | None): pressure - biotCoefficient (npt.NDArray[np.float64]): biot coefficient - totalStressAttributeName (str): total stress attribute name - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - totalStress: npt.NDArray[ np.float64 ] - # if pressure data is missing, total stress equals effective stress - if pressure is None: - # copy effective stress data - totalStress = np.copy( effectiveStress ) - mess: str = ( "Pressure attribute is undefined, " + "total stress will be equal to effective stress." ) - self.m_logger.warning( mess ) - else: - if np.isnan( pressure ).any(): - mess0: str = ( "Some cells do not have pressure data, " + - "total stress will be equal to effective stress." ) - self.m_logger.warning( mess0 ) - - # replace nan values by 0. - pressure[ np.isnan( pressure ) ] = 0.0 - totalStress = fcts.totalStress( effectiveStress, biotCoefficient, pressure ) - - # create total stress attribute - assert totalStress is not None, "Total stress data is null." - createAttribute( - self.m_output, - totalStress, - totalStressAttributeName, - ComponentNameEnum.XYZ.value, - self.m_attributeOnPoints, - ) - return True - - def computeLitostaticStress( self: Self ) -> bool: - """Compute lithostatic stress. - - Returns: - bool: True if calculation successfully ended, False otherwise. - """ - litostaticStressAttributeName: str = ( PostProcessingOutputsEnum.LITHOSTATIC_STRESS.attributeName ) - if not isAttributeInObject( self.m_output, litostaticStressAttributeName, self.m_attributeOnPoints ): - - densityAttributeName: str = GeosMeshOutputsEnum.ROCK_DENSITY.attributeName - density: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, densityAttributeName, - self.m_attributeOnPoints ) - try: - depth: npt.NDArray[ np.float64 ] = self.computeDepthAlongLine( - ) if self.m_attributeOnPoints else self.computeDepthInMesh() - assert depth is not None, "Depth is undefined." - assert density is not None, ( f"{densityAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - litostaticStress = fcts.lithostaticStress( depth, density, GRAVITY ) - createAttribute( - self.m_output, - litostaticStress, - litostaticStressAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Lithostatic stress was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeDepthAlongLine( self: Self ) -> npt.NDArray[ np.float64 ]: - """Compute depth along a line. - - Returns: - npt.NDArray[np.float64]: 1D array with depth property - """ - # get z coordinate - zCoord: npt.NDArray[ np.float64 ] = self.getZcoordinates() - assert zCoord is not None, "Depth coordinates cannot be computed." - - # TODO: to find how to compute depth in a general case - # GEOS z axis is upward - depth: npt.NDArray[ np.float64 ] = -1.0 * zCoord - return depth - - def computeDepthInMesh( self: Self ) -> npt.NDArray[ np.float64 ]: - """Compute depth of each cell in a mesh. - - Returns: - npt.NDArray[np.float64]: array with depth property - """ - # get z coordinate - zCoord: npt.NDArray[ np.float64 ] = self.getZcoordinates() - assert zCoord is not None, "Depth coordinates cannot be computed." - - # TODO: to find how to compute depth in a general case - depth: npt.NDArray[ np.float64 ] = -1.0 * zCoord - return depth - - def getZcoordinates( self: Self ) -> npt.NDArray[ np.float64 ]: - """Get z coordinates from self.m_output. - - Returns: - npt.NDArray[np.float64]: 1D array with depth property - """ - # get z coordinate - zCoord: npt.NDArray[ np.float64 ] - pointCoords: npt.NDArray[ np.float64 ] = self.getPointCoordinates() - assert pointCoords is not None, "Point coordinates are undefined." - assert pointCoords.shape[ 1 ] == 2, "Point coordinates are undefined." - zCoord = pointCoords[ :, 2 ] - return zCoord - - def computeElasticStrain( self: Self ) -> bool: - """Compute elastic strain from effective stress and elastic modulus. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - elasticStrainAttributeName: str = ( PostProcessingOutputsEnum.STRAIN_ELASTIC.attributeName ) - if not isAttributeInObject( self.m_output, elasticStrainAttributeName, self.m_attributeOnPoints ): - effectiveStressAttributeName: str = ( GeosMeshOutputsEnum.STRESS_EFFECTIVE.attributeName ) - effectiveStress: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, effectiveStressAttributeName, - self.m_attributeOnPoints ) - - effectiveStressIniAttributeName: str = ( PostProcessingOutputsEnum.STRESS_EFFECTIVE_INITIAL.attributeName ) - effectiveStressIni: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, - effectiveStressIniAttributeName, - self.m_attributeOnPoints ) - - bulkModulusAttributeName: str = ( GeosMeshOutputsEnum.BULK_MODULUS.attributeName ) - bulkModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, bulkModulusAttributeName, - self.m_attributeOnPoints ) - - shearModulusAttributeName: str = ( GeosMeshOutputsEnum.SHEAR_MODULUS.attributeName ) - shearModulus: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, shearModulusAttributeName, - self.m_attributeOnPoints ) - - try: - assert effectiveStress is not None, ( f"{effectiveStressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert effectiveStressIni is not None, ( f"{effectiveStressIniAttributeName}" + - UNDEFINED_ATTRIBUTE_MESSAGE ) - assert bulkModulus is not None, ( f"{bulkModulusAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert shearModulus is not None, ( f"{shearModulusAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - deltaEffectiveStress = effectiveStress - effectiveStressIni - elasticStrain: npt.NDArray[ np.float64 ] = ( fcts.elasticStrainFromBulkShear( - deltaEffectiveStress, bulkModulus, shearModulus ) ) - createAttribute( - self.m_output, - elasticStrain, - elasticStrainAttributeName, - ComponentNameEnum.XYZ.value, - self.m_attributeOnPoints, - ) - - except AssertionError as e: - self.m_logger.error( "Elastic strain was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - return True - - def computeReservoirStressPathReal( self: Self ) -> bool: - """Compute reservoir stress paths. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - RSPrealAttributeName: str = PostProcessingOutputsEnum.RSP_REAL.attributeName - if not isAttributeInObject( self.m_output, RSPrealAttributeName, self.m_attributeOnPoints ): - # real RSP - try: - # total stress at current and initial time steps - totalStressAttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL.attributeName ) - totalStress: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, totalStressAttributeName, - self.m_attributeOnPoints ) - - totalStressT0AttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL_INITIAL.attributeName ) - totalStressIni: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, totalStressT0AttributeName, - self.m_attributeOnPoints ) - - deltaPressureAttributeName: str = ( GeosMeshOutputsEnum.DELTA_PRESSURE.attributeName ) - deltaPressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, deltaPressureAttributeName, - self.m_attributeOnPoints ) - - assert totalStress is not None, ( f"{totalStressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert totalStressIni is not None, ( f"{totalStressT0AttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert deltaPressure is not None, ( f"{deltaPressureAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - # create delta stress attribute for QC - deltaTotalStressAttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL_DELTA.attributeName ) - deltaStress: npt.NDArray[ np.float64 ] = totalStress - totalStressIni - componentNames: tuple[ str, ...] = getComponentNames( self.m_output, totalStressAttributeName, - self.m_attributeOnPoints ) - assert deltaStress is not None, ( f"{deltaTotalStressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - createAttribute( - self.m_output, - deltaStress, - deltaTotalStressAttributeName, - componentNames, - self.m_attributeOnPoints, - ) - - RSPreal: npt.NDArray[ np.float64 ] = fcts.reservoirStressPathReal( deltaStress, deltaPressure ) - assert RSPreal is not None, ( f"{RSPrealAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - createAttribute( - self.m_output, - RSPreal, - RSPrealAttributeName, - componentNames, - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Reservoir stress path in real conditions was " + "not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeReservoirStressPathOed( self: Self ) -> bool: - """Compute Reservoir Stress Path in oedometric conditions. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - RSPOedAttributeName: str = PostProcessingOutputsEnum.RSP_OED.attributeName - if not isAttributeInObject( self.m_output, RSPOedAttributeName, self.m_attributeOnPoints ): - - poissonRatioAttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO.attributeName ) - poissonRatio: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, poissonRatioAttributeName, - self.m_attributeOnPoints ) - - biotCoefficientAttributeName: str = ( PostProcessingOutputsEnum.BIOT_COEFFICIENT.attributeName ) - biotCoefficient: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, biotCoefficientAttributeName, - self.m_attributeOnPoints ) - - try: - assert poissonRatio is not None, ( f"{poissonRatioAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert biotCoefficient is not None, ( f"{biotCoefficientAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - RSPoed: npt.NDArray[ np.float64 ] = fcts.reservoirStressPathOed( biotCoefficient, poissonRatio ) - createAttribute( - self.m_output, - RSPoed, - RSPOedAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Oedometric RSP was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeStressRatioReal( self: Self, inputAttribute: AttributeEnum, outputAttribute: AttributeEnum ) -> bool: - """Compute the ratio between horizontal and vertical effective stress. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - stressRatioRealAttributeName = outputAttribute.attributeName - if not isAttributeInObject( self.m_output, stressRatioRealAttributeName, self.m_attributeOnPoints ): - try: - stressAttributeName: str = inputAttribute.attributeName - stress: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, stressAttributeName, - self.m_attributeOnPoints ) - assert stress is not None, ( f"{stressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - verticalStress: npt.NDArray[ np.float64 ] = stress[ :, 2 ] - # keep the minimum of the 2 horizontal components - horizontalStress: npt.NDArray[ np.float64 ] = np.min( stress[ :, :2 ], axis=1 ) - - stressRatioReal: npt.NDArray[ np.float64 ] = fcts.stressRatio( horizontalStress, verticalStress ) - createAttribute( - self.m_output, - stressRatioReal, - stressRatioRealAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( f"{outputAttribute.attributeName} was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - return True - - def computeEffectiveStressRatioOed( self: Self ) -> bool: - """Compute the effective stress ratio in oedometric conditions. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - effectiveStressRatioOedAttributeName: str = ( - PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_OED.attributeName ) - if not isAttributeInObject( - self.m_output, - effectiveStressRatioOedAttributeName, - self.m_attributeOnPoints, - ): - poissonRatioAttributeName: str = ( PostProcessingOutputsEnum.POISSON_RATIO.attributeName ) - poissonRatio: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, poissonRatioAttributeName, - self.m_attributeOnPoints ) - - try: - assert poissonRatio is not None, ( f"{poissonRatioAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - - effectiveStressRatioOed: npt.NDArray[ np.float64 ] = ( fcts.deviatoricStressPathOed( poissonRatio ) ) - createAttribute( - self.m_output, - effectiveStressRatioOed, - effectiveStressRatioOedAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Oedometric effective stress ratio was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - return True - - def computeCriticalTotalStressRatio( self: Self ) -> bool: - """Compute fracture index and fracture threshold. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - ret: bool = True - - fractureIndexAttributeName = ( PostProcessingOutputsEnum.CRITICAL_TOTAL_STRESS_RATIO.attributeName ) - if not isAttributeInObject( self.m_output, fractureIndexAttributeName, self.m_attributeOnPoints ): - - stressAttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL.attributeName ) - stress: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, stressAttributeName, - self.m_attributeOnPoints ) - - pressureAttributeName: str = GeosMeshOutputsEnum.PRESSURE.attributeName - pressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, pressureAttributeName, - self.m_attributeOnPoints ) - - try: - assert stress is not None, ( f"{stressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert pressure is not None, ( f"{pressureAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert stress.shape[ 1 ] >= 3 - except AssertionError as e: - self.m_logger.error( "Critical total stress ratio and threshold were not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - # fracture index - try: - verticalStress: npt.NDArray[ np.float64 ] = stress[ :, 2 ] - criticalTotalStressRatio: npt.NDArray[ np.float64 ] = ( fcts.criticalTotalStressRatio( - pressure, verticalStress ) ) - createAttribute( - self.m_output, - criticalTotalStressRatio, - fractureIndexAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Critical total stress ratio was not computed due to:" ) - self.m_logger.error( str( e ) ) - ret = False - - # fracture threshold - try: - mask: npt.NDArray[ np.bool_ ] = np.argmin( np.abs( stress[ :, :2 ] ), axis=1 ) - horizontalStress: npt.NDArray[ np.float64 ] = stress[ :, :2 ][ np.arange( stress[ :, :2 ].shape[ 0 ] ), - mask ] - - stressRatioThreshold: npt.NDArray[ np.float64 ] = ( fcts.totalStressRatioThreshold( - pressure, horizontalStress ) ) - fractureThresholdAttributeName = ( - PostProcessingOutputsEnum.TOTAL_STRESS_RATIO_THRESHOLD.attributeName ) - createAttribute( - self.m_output, - stressRatioThreshold, - fractureThresholdAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Total stress ratio threshold was not computed due to:" ) - self.m_logger.error( str( e ) ) - ret = False - - return ret - - def computeCriticalPorePressure( self: Self ) -> bool: - """Compute the critical pore pressure and the pressure index. - - Returns: - bool: return True if calculation successfully ended, False otherwise. - """ - ret: bool = True - criticalPorePressureAttributeName = ( PostProcessingOutputsEnum.CRITICAL_PORE_PRESSURE.attributeName ) - if not isAttributeInObject( self.m_output, criticalPorePressureAttributeName, self.m_attributeOnPoints ): - - stressAttributeName: str = ( PostProcessingOutputsEnum.STRESS_TOTAL.attributeName ) - stress: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, stressAttributeName, - self.m_attributeOnPoints ) - - pressureAttributeName: str = GeosMeshOutputsEnum.PRESSURE.attributeName - pressure: npt.NDArray[ np.float64 ] = getArrayInObject( self.m_output, pressureAttributeName, - self.m_attributeOnPoints ) - - try: - assert self.m_rockCohesion is not None, "Rock cohesion is undefined." - assert self.m_frictionAngle is not None, "Friction angle is undefined." - assert stress is not None, ( f"{stressAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert pressure is not None, ( f"{pressureAttributeName}" + UNDEFINED_ATTRIBUTE_MESSAGE ) - assert stress.shape[ 1 ] >= 3 - except AssertionError as e: - self.m_logger.error( "Critical pore pressure and threshold were not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - try: - criticalPorePressure: npt.NDArray[ np.float64 ] = ( fcts.criticalPorePressure( - -1.0 * stress, self.m_rockCohesion, self.m_frictionAngle ) ) - createAttribute( - self.m_output, - criticalPorePressure, - criticalPorePressureAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Critical pore pressure was not computed due to:" ) - self.m_logger.error( str( e ) ) - return False - - try: - # add critical pore pressure index (i.e., ratio between pressure and criticalPorePressure) - assert criticalPorePressure is not None, ( "Critical pore pressure attribute" + " was not created" ) - criticalPorePressureIndex: npt.NDArray[ np.float64 ] = ( fcts.criticalPorePressureThreshold( - pressure, criticalPorePressure ) ) - criticalPorePressureIndexAttributeName = ( - PostProcessingOutputsEnum.CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName ) - createAttribute( - self.m_output, - criticalPorePressureIndex, - criticalPorePressureIndexAttributeName, - (), - self.m_attributeOnPoints, - ) - except AssertionError as e: - self.m_logger.error( "Pore pressure threshold was not computed due to:" ) - self.m_logger.error( str( e ) ) - ret = False - - return ret - - def getPointCoordinates( self: Self ) -> npt.NDArray[ np.float64 ]: - """Get the coordinates of Points or Cell center. - - Returns: - npt.NDArray[np.float64]: points/cell center coordinates - """ - if self.m_attributeOnPoints: - return self.m_output.GetPoints() # type: ignore[no-any-return] - else: - # Find cell centers - filter = vtkCellCenters() - filter.SetInputDataObject( self.m_output ) - filter.Update() - return filter.GetOutput().GetPoints() # type: ignore[no-any-return] diff --git a/geos-processing/pyproject.toml b/geos-processing/pyproject.toml new file mode 100644 index 00000000..71fadf3e --- /dev/null +++ b/geos-processing/pyproject.toml @@ -0,0 +1,70 @@ +[build-system] +requires = ["setuptools>=61.2", "wheel >= 0.37.1"] +build-backend = "setuptools.build_meta" + +[tool.setuptools] +include-package-data = true + +[tool.setuptools.packages.find] +where = ["src"] +include = ["geos.processing*"] +exclude = ['tests*'] + +[project] +name = "geos-processing" +version = "0.1.0" +description = "GEOS post- and pre-processing tools" +authors = [{name = "GEOS Contributors" }] +maintainers = [ + {name = "Alexandre Benedicto", email = "alexandre.benedicto@external.totalenergies.com" }, + {name = "Romain Baville", email = "romain.baville@external.totalenergies.com" }, + {name = "Paloma Martinez", email = "paloma.martinez@external.totalenergies.com" }, +] +license = {text = "Apache-2.0"} +classifiers = [ + "Development Status :: 4 - Beta", + "Programming Language :: Python" +] + +requires-python = ">=3.10" + +dependencies = [ + "geos-geomechanics", + "geos-utils", + "geos-mesh", + "vtk >= 9.3", + "numpy >= 2.2", + "typing_extensions >= 4.12", +] + + +[project.urls] +Homepage = "https://github.com/GEOS-DEV/geosPythonPackages" +Documentation = "https://geosx-geosx.readthedocs-hosted.com/projects/geosx-geospythonpackages/en/latest/" +Repository = "https://github.com/GEOS-DEV/geosPythonPackages.git" +"Bug Tracker" = "https://github.com/GEOS-DEV/geosPythonPackages/issues" + +[project.optional-dependencies] +build = [ + "build ~= 1.2" +] +dev = [ + "mypy", + "ruff", + "yapf", +] +test = [ + "pytest-cov", + "pytest" +] + +[tool.pytest.ini_options] +addopts = "--import-mode=importlib" +console_output_style = "count" +pythonpath = ["src"] +python_classes = "Test" +python_files = "test*.py" +python_functions = "test*" +testpaths = ["tests"] +norecursedirs = "bin" +filterwarnings = [] \ No newline at end of file diff --git a/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py b/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py new file mode 100644 index 00000000..ecef2cfc --- /dev/null +++ b/geos-processing/src/geos/processing/post_processing/GeomechanicsCalculator.py @@ -0,0 +1,1454 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Martin Lemay, Romain Baville +# ruff: noqa: E402 # disable Module level import not at top of file +from typing import Union +from typing_extensions import Self +from dataclasses import dataclass + +import logging +import numpy as np +import numpy.typing as npt + +import geos.geomechanics.processing.geomechanicsCalculatorFunctions as fcts + +from geos.mesh.utils.arrayModifiers import createAttribute +from geos.mesh.utils.arrayHelpers import ( + getArrayInObject, + isAttributeInObject, +) + +from geos.utils.Logger import ( + Logger, + getLogger, +) +from geos.utils.GeosOutputsConstants import ( + AttributeEnum, + ComponentNameEnum, + GeosMeshOutputsEnum, + PostProcessingOutputsEnum, +) +from geos.utils.PhysicalConstants import ( + DEFAULT_FRICTION_ANGLE_RAD, + DEFAULT_GRAIN_BULK_MODULUS, + DEFAULT_ROCK_COHESION, + GRAVITY, + WATER_DENSITY, +) + +from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid +from vtkmodules.vtkFiltersCore import vtkCellCenters + +__doc__ = """ +GeomechanicsCalculator module is a vtk filter that allows to compute additional geomechanical properties from existing ones: + - The Young modulus and the Poisson's ratio named "youngModulus" and "poissonRatio" or bulk and shear moduli named "bulkModulus" and "shearModulus" + - The initial Young modulus and Poisson's ratio named "youngModulusInitial" and "poissonRatioInitial" or the initial bulk modulus named "bulkModulusInitial" + - The porosity named "porosity" + - The initial porosity named "porosityInitial" + - The delta of pressure named "deltaPressure" + - The density named "density" + - The effective stress named "stressEffective" + - The initial effective stress named "stressEffectiveInitial" + - The pressure named "pressure" + +The basic geomechanics outputs are: + - The elastic moduli not present on the mesh + - Biot coefficient + - Compressibility, oedometric compressibility and real compressibility coefficient + - Specific gravity + - Real effective stress ratio + - Total initial stress, total current stress and total stress ratio + - Lithostatic stress (physic to update) + - Elastic stain + - Real reservoir stress path and reservoir stress path in oedometric condition + +The advanced geomechanics outputs are: + - Fracture index and threshold + - Critical pore pressure and pressure index + +GeomechanicsCalculator filter input mesh is either vtkPointSet or vtkUnstructuredGrid +and returned mesh is of same type as input. + +.. Important:: + Please refer to the GeosExtractMergeBlockVolume* filters to provide the correct input. + +.. Note:: + - The default physical constants used by the filter are: + - grainBulkModulus = 38e9 Pa ( quartz value ) + - specificDensity = 1000.0 kg/m³ ( water value ) + - rockCohesion = 0.0 Pa ( fractured case ) + - frictionAngle = 10.0° + +To use the filter: + +.. code-block:: python + + import numpy as np + from geos.mesh.processing.GeomechanicsCalculator import GeomechanicsCalculator + + # Define filter inputs + mesh: Union[ vtkPointSet, vtkUnstructuredGrid ] + computeAdvancedOutputs: bool # optional, defaults to False + speHandler: bool # optional, defaults to False + + # Instantiate the filter + filter: GeomechanicsCalculator = GeomechanicsCalculator( mesh, computeAdvancedOutputs, speHandler ) + + # Use your own handler (if speHandler is True) + yourHandler: logging.Handler + filter.setLoggerHandler( yourHandler ) + + # Change the physical constants if needed + ## Basic outputs + grainBulkModulus: float + specificDensity: float + filter.physicalConstants.grainBulkModulus = grainBulkModulus + filter.physicalConstants.specificDensity = specificDensity + + ## Advanced outputs + rockCohesion: float + frictionAngle: float + filter.physicalConstants.rockCohesion = rockCohesion + filter.physicalConstants.frictionAngle = frictionAngle + + # Do calculations + filter.applyFilter() + + # Get the mesh with the geomechanical output as attribute + output: Union[vtkPointSet, vtkUnstructuredGrid] + output = filter.getOutput() +""" + +loggerTitle: str = "Geomechanical Calculator Filter" + +# Elastic Moduli: +BULK_MODULUS: AttributeEnum = GeosMeshOutputsEnum.BULK_MODULUS +SHEAR_MODULUS: AttributeEnum = GeosMeshOutputsEnum.SHEAR_MODULUS +YOUNG_MODULUS: AttributeEnum = PostProcessingOutputsEnum.YOUNG_MODULUS +POISSON_RATIO: AttributeEnum = PostProcessingOutputsEnum.POISSON_RATIO +BULK_MODULUS_T0: AttributeEnum = PostProcessingOutputsEnum.BULK_MODULUS_INITIAL +YOUNG_MODULUS_T0: AttributeEnum = PostProcessingOutputsEnum.YOUNG_MODULUS_INITIAL +POISSON_RATIO_T0: AttributeEnum = PostProcessingOutputsEnum.POISSON_RATIO_INITIAL +ELASTIC_MODULI: tuple[ AttributeEnum, ...] = ( BULK_MODULUS, SHEAR_MODULUS, YOUNG_MODULUS, POISSON_RATIO, + BULK_MODULUS_T0, YOUNG_MODULUS_T0, POISSON_RATIO_T0 ) + +# Mandatory attributes: +POROSITY: AttributeEnum = GeosMeshOutputsEnum.POROSITY +POROSITY_T0: AttributeEnum = GeosMeshOutputsEnum.POROSITY_INI +PRESSURE: AttributeEnum = GeosMeshOutputsEnum.PRESSURE +DELTA_PRESSURE: AttributeEnum = GeosMeshOutputsEnum.DELTA_PRESSURE +DENSITY: AttributeEnum = GeosMeshOutputsEnum.ROCK_DENSITY +STRESS_EFFECTIVE: AttributeEnum = GeosMeshOutputsEnum.STRESS_EFFECTIVE +STRESS_EFFECTIVE_T0: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_INITIAL +MANDATORY_ATTRIBUTES: tuple[ AttributeEnum, ...] = ( POROSITY, POROSITY_T0, PRESSURE, DELTA_PRESSURE, DENSITY, + STRESS_EFFECTIVE, STRESS_EFFECTIVE_T0 ) + +# Basic outputs: +BIOT_COEFFICIENT: AttributeEnum = PostProcessingOutputsEnum.BIOT_COEFFICIENT +COMPRESSIBILITY: AttributeEnum = PostProcessingOutputsEnum.COMPRESSIBILITY +COMPRESSIBILITY_OED: AttributeEnum = PostProcessingOutputsEnum.COMPRESSIBILITY_OED +COMPRESSIBILITY_REAL: AttributeEnum = PostProcessingOutputsEnum.COMPRESSIBILITY_REAL +SPECIFIC_GRAVITY: AttributeEnum = PostProcessingOutputsEnum.SPECIFIC_GRAVITY +STRESS_EFFECTIVE_RATIO_REAL: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_REAL +STRESS_TOTAL: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL +STRESS_TOTAL_T0: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL_INITIAL +STRESS_TOTAL_RATIO_REAL: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL_RATIO_REAL +LITHOSTATIC_STRESS: AttributeEnum = PostProcessingOutputsEnum.LITHOSTATIC_STRESS +STRAIN_ELASTIC: AttributeEnum = PostProcessingOutputsEnum.STRAIN_ELASTIC +STRESS_TOTAL_DELTA: AttributeEnum = PostProcessingOutputsEnum.STRESS_TOTAL_DELTA +RSP_REAL: AttributeEnum = PostProcessingOutputsEnum.RSP_REAL +RSP_OED: AttributeEnum = PostProcessingOutputsEnum.RSP_OED +STRESS_EFFECTIVE_RATIO_OED: AttributeEnum = PostProcessingOutputsEnum.STRESS_EFFECTIVE_RATIO_OED +BASIC_OUTPUTS: tuple[ AttributeEnum, + ...] = ( BIOT_COEFFICIENT, COMPRESSIBILITY, COMPRESSIBILITY_OED, COMPRESSIBILITY_REAL, + SPECIFIC_GRAVITY, STRESS_EFFECTIVE_RATIO_REAL, STRESS_TOTAL, STRESS_TOTAL_T0, + STRESS_TOTAL_RATIO_REAL, LITHOSTATIC_STRESS, STRAIN_ELASTIC, STRESS_TOTAL_DELTA, + RSP_REAL, RSP_OED, STRESS_EFFECTIVE_RATIO_OED ) + +# Advanced outputs: +CRITICAL_TOTAL_STRESS_RATIO: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_TOTAL_STRESS_RATIO +TOTAL_STRESS_RATIO_THRESHOLD: AttributeEnum = PostProcessingOutputsEnum.TOTAL_STRESS_RATIO_THRESHOLD +CRITICAL_PORE_PRESSURE: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_PORE_PRESSURE +CRITICAL_PORE_PRESSURE_THRESHOLD: AttributeEnum = PostProcessingOutputsEnum.CRITICAL_PORE_PRESSURE_THRESHOLD +ADVANCED_OUTPUTS: tuple[ AttributeEnum, ...] = ( CRITICAL_TOTAL_STRESS_RATIO, TOTAL_STRESS_RATIO_THRESHOLD, + CRITICAL_PORE_PRESSURE, CRITICAL_PORE_PRESSURE_THRESHOLD ) + + +class GeomechanicsCalculator: + + @dataclass + class PhysicalConstants: + """The dataclass with the value of the physical constant used to compute geomechanics properties.""" + ## Basic outputs + _grainBulkModulus: float = DEFAULT_GRAIN_BULK_MODULUS + _specificDensity: float = WATER_DENSITY + ## Advanced outputs + _rockCohesion: float = DEFAULT_ROCK_COHESION + _frictionAngle: float = DEFAULT_FRICTION_ANGLE_RAD + + @property + def grainBulkModulus( self: Self ) -> float: + """Get the grain bulk modulus value.""" + return self._grainBulkModulus + + @grainBulkModulus.setter + def grainBulkModulus( self: Self, value: float ) -> None: + self._grainBulkModulus = value + + @property + def specificDensity( self: Self ) -> float: + """Get the specific density value.""" + return self._specificDensity + + @specificDensity.setter + def specificDensity( self: Self, value: float ) -> None: + self._specificDensity = value + + @property + def rockCohesion( self: Self ) -> float: + """Get the rock cohesion value.""" + return self._rockCohesion + + @rockCohesion.setter + def rockCohesion( self: Self, value: float ) -> None: + self._rockCohesion = value + + @property + def frictionAngle( self: Self ) -> float: + """Get the friction angle value.""" + return self._frictionAngle + + @frictionAngle.setter + def frictionAngle( self: Self, value: float ) -> None: + self._frictionAngle = value + + @dataclass + class ElasticModuliValue: + """The dataclass with the value of the elastic moduli.""" + _bulkModulus: npt.NDArray[ np.float64 ] | None = None + _shearModulus: npt.NDArray[ np.float64 ] | None = None + _youngModulus: npt.NDArray[ np.float64 ] | None = None + _poissonRatio: npt.NDArray[ np.float64 ] | None = None + _bulkModulusT0: npt.NDArray[ np.float64 ] | None = None + _youngModulusT0: npt.NDArray[ np.float64 ] | None = None + _poissonRatioT0: npt.NDArray[ np.float64 ] | None = None + + @property + def bulkModulus( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the bulk modulus value.""" + return self._bulkModulus + + @bulkModulus.setter + def bulkModulus( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._bulkModulus = value + + @property + def shearModulus( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the shear modulus value.""" + return self._shearModulus + + @shearModulus.setter + def shearModulus( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._shearModulus = value + + @property + def youngModulus( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the young modulus value.""" + return self._youngModulus + + @youngModulus.setter + def youngModulus( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._youngModulus = value + + @property + def poissonRatio( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the poisson ratio value.""" + return self._poissonRatio + + @poissonRatio.setter + def poissonRatio( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._poissonRatio = value + + @property + def bulkModulusT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the bulk modulus at the initial time value.""" + return self._bulkModulusT0 + + @bulkModulusT0.setter + def bulkModulusT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._bulkModulusT0 = value + + @property + def youngModulusT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the young modulus at the initial time value.""" + return self._youngModulusT0 + + @youngModulusT0.setter + def youngModulusT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._youngModulusT0 = value + + @property + def poissonRatioT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the poisson ration at the initial time value.""" + return self._poissonRatioT0 + + @poissonRatioT0.setter + def poissonRatioT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._poissonRatioT0 = value + + def setElasticModulusValue( self: Self, name: str, value: npt.NDArray[ np.float64 ] ) -> None: + """Set the elastic modulus value wanted. + + Args: + name (str): The name of the elastic modulus. + value (npt.NDArray[np.float64]): The value to set. + """ + if name == BULK_MODULUS.attributeName: + self.bulkModulus = value + elif name == BULK_MODULUS_T0.attributeName: + self.bulkModulusT0 = value + elif name == SHEAR_MODULUS.attributeName: + self.shearModulus = value + elif name == YOUNG_MODULUS.attributeName: + self.youngModulus = value + elif name == YOUNG_MODULUS_T0.attributeName: + self.youngModulusT0 = value + elif name == POISSON_RATIO.attributeName: + self.poissonRatio = value + elif name == POISSON_RATIO_T0.attributeName: + self.poissonRatioT0 = value + + def getElasticModulusValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ] | None: + """Get the wanted elastic modulus value. + + Args: + name (str): The name of the wanted elastic modulus. + + Returns: + npt.NDArray[np.float64]: The value of the elastic modulus. + """ + if name == BULK_MODULUS.attributeName: + return self.bulkModulus + elif name == BULK_MODULUS_T0.attributeName: + return self.bulkModulusT0 + elif name == SHEAR_MODULUS.attributeName: + return self.shearModulus + elif name == YOUNG_MODULUS.attributeName: + return self.youngModulus + elif name == YOUNG_MODULUS_T0.attributeName: + return self.youngModulusT0 + elif name == POISSON_RATIO.attributeName: + return self.poissonRatio + elif name == POISSON_RATIO_T0.attributeName: + return self.poissonRatioT0 + else: + raise NameError + + @dataclass + class MandatoryAttributesValue: + """The dataclass with the value of mandatory properties to have to compute other geomechanics properties.""" + _porosity: npt.NDArray[ np.float64 ] | None = None + _porosityInitial: npt.NDArray[ np.float64 ] | None = None + _pressure: npt.NDArray[ np.float64 ] | None = None + _deltaPressure: npt.NDArray[ np.float64 ] | None = None + _density: npt.NDArray[ np.float64 ] | None = None + _effectiveStress: npt.NDArray[ np.float64 ] | None = None + _effectiveStressT0: npt.NDArray[ np.float64 ] | None = None + + @property + def porosity( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the porosity value.""" + return self._porosity + + @property + def porosityInitial( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the initial porosity value.""" + return self._porosityInitial + + @property + def pressure( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the pressure value.""" + return self._pressure + + @property + def deltaPressure( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the delta pressure value.""" + return self._deltaPressure + + @property + def density( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the density value.""" + return self._density + + @property + def effectiveStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the effective stress value.""" + return self._effectiveStress + + @property + def effectiveStressT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the initial effective stress value.""" + return self._effectiveStressT0 + + def setMandatoryAttributeValue( self: Self, name: str, value: npt.NDArray[ np.float64 ] ) -> None: + """Set the value of a mandatory property. + + Args: + name (str): The name of the property. + value (npt.NDArray[np.float64]): The value to set. + """ + if name == POROSITY.attributeName: + self._porosity = value + elif name == POROSITY_T0.attributeName: + self._porosityInitial = value + elif name == PRESSURE.attributeName: + self._pressure = value + elif name == DELTA_PRESSURE.attributeName: + self._deltaPressure = value + elif name == DENSITY.attributeName: + self._density = value + elif name == STRESS_EFFECTIVE.attributeName: + self._effectiveStress = value + elif name == STRESS_EFFECTIVE_T0.attributeName: + self._effectiveStressT0 = value + + @dataclass + class BasicOutputValue: + """The dataclass with the value of the basic geomechanics outputs.""" + _biotCoefficient: npt.NDArray[ np.float64 ] | None = None + _compressibility: npt.NDArray[ np.float64 ] | None = None + _compressibilityOed: npt.NDArray[ np.float64 ] | None = None + _compressibilityReal: npt.NDArray[ np.float64 ] | None = None + _specificGravity: npt.NDArray[ np.float64 ] | None = None + _effectiveStressRatioReal: npt.NDArray[ np.float64 ] | None = None + _totalStress: npt.NDArray[ np.float64 ] | None = None + _totalStressT0: npt.NDArray[ np.float64 ] | None = None + _totalStressRatioReal: npt.NDArray[ np.float64 ] | None = None + # _lithostaticStress: npt.NDArray[ np.float64 ] | None = None + _elasticStrain: npt.NDArray[ np.float64 ] | None = None + _deltaTotalStress: npt.NDArray[ np.float64 ] | None = None + _rspReal: npt.NDArray[ np.float64 ] | None = None + _rspOed: npt.NDArray[ np.float64 ] | None = None + _effectiveStressRatioOed: npt.NDArray[ np.float64 ] | None = None + + @property + def biotCoefficient( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the biot coefficient value.""" + return self._biotCoefficient + + @biotCoefficient.setter + def biotCoefficient( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._biotCoefficient = value + + @property + def compressibility( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the compressibility value.""" + return self._compressibility + + @compressibility.setter + def compressibility( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._compressibility = value + + @property + def compressibilityOed( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the compressibility in oedometric condition value.""" + return self._compressibilityOed + + @compressibilityOed.setter + def compressibilityOed( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._compressibilityOed = value + + @property + def compressibilityReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the real compressibility value.""" + return self._compressibilityReal + + @compressibilityReal.setter + def compressibilityReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._compressibilityReal = value + + @property + def specificGravity( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the specific gravity value.""" + return self._specificGravity + + @specificGravity.setter + def specificGravity( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._specificGravity = value + + @property + def effectiveStressRatioReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the real effective stress ratio value.""" + return self._effectiveStressRatioReal + + @effectiveStressRatioReal.setter + def effectiveStressRatioReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._effectiveStressRatioReal = value + + @property + def totalStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the total stress value.""" + return self._totalStress + + @totalStress.setter + def totalStress( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._totalStress = value + + @property + def totalStressT0( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the initial total stress value.""" + return self._totalStressT0 + + @totalStressT0.setter + def totalStressT0( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._totalStressT0 = value + + @property + def totalStressRatioReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the total real stress ratio value.""" + return self._totalStressRatioReal + + @totalStressRatioReal.setter + def totalStressRatioReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._totalStressRatioReal = value + + @property + def lithostaticStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the lithostatic stress value.""" + return self._lithostaticStress + + @lithostaticStress.setter + def lithostaticStress( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._lithostaticStress = value + + @property + def elasticStrain( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the elastic strain value.""" + return self._elasticStrain + + @elasticStrain.setter + def elasticStrain( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._elasticStrain = value + + @property + def deltaTotalStress( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the total delta stress value.""" + return self._deltaTotalStress + + @deltaTotalStress.setter + def deltaTotalStress( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._deltaTotalStress = value + + @property + def rspReal( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the real reservoir stress path value.""" + return self._rspReal + + @rspReal.setter + def rspReal( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._rspReal = value + + @property + def rspOed( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the reservoir stress path in oedometric condition value.""" + return self._rspOed + + @rspOed.setter + def rspOed( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._rspOed = value + + @property + def effectiveStressRatioOed( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the effective stress ratio oedometric value.""" + return self._effectiveStressRatioOed + + @effectiveStressRatioOed.setter + def effectiveStressRatioOed( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._effectiveStressRatioOed = value + + def getBasicOutputValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ] | None: + """Get the value of the basic output wanted. + + Args: + name (str): The name of the basic output. + + Returns: + npt.NDArray[ np.float64 ]: the value of the basic output. + """ + if name == BIOT_COEFFICIENT.attributeName: + return self.biotCoefficient + elif name == COMPRESSIBILITY.attributeName: + return self.compressibility + elif name == COMPRESSIBILITY_OED.attributeName: + return self.compressibilityOed + elif name == COMPRESSIBILITY_REAL.attributeName: + return self.compressibilityReal + elif name == SPECIFIC_GRAVITY.attributeName: + return self.specificGravity + elif name == STRESS_EFFECTIVE_RATIO_REAL.attributeName: + return self.effectiveStressRatioReal + elif name == STRESS_TOTAL.attributeName: + return self.totalStress + elif name == STRESS_TOTAL_T0.attributeName: + return self.totalStressT0 + elif name == STRESS_TOTAL_RATIO_REAL.attributeName: + return self.totalStressRatioReal + elif name == LITHOSTATIC_STRESS.attributeName: + return self.lithostaticStress + elif name == STRAIN_ELASTIC.attributeName: + return self.elasticStrain + elif name == STRESS_TOTAL_DELTA.attributeName: + return self.deltaTotalStress + elif name == RSP_REAL.attributeName: + return self.rspReal + elif name == RSP_OED.attributeName: + return self.rspOed + elif name == STRESS_EFFECTIVE_RATIO_OED.attributeName: + return self.effectiveStressRatioOed + else: + raise NameError + + @dataclass + class AdvancedOutputValue: + """The dataclass with the value of the advanced geomechanics outputs.""" + _criticalTotalStressRatio: npt.NDArray[ np.float64 ] | None = None + _stressRatioThreshold: npt.NDArray[ np.float64 ] | None = None + _criticalPorePressure: npt.NDArray[ np.float64 ] | None = None + _criticalPorePressureIndex: npt.NDArray[ np.float64 ] | None = None + + @property + def criticalTotalStressRatio( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the critical total stress ratio value.""" + return self._criticalTotalStressRatio + + @criticalTotalStressRatio.setter + def criticalTotalStressRatio( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._criticalTotalStressRatio = value + + @property + def stressRatioThreshold( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the stress ratio threshold value.""" + return self._stressRatioThreshold + + @stressRatioThreshold.setter + def stressRatioThreshold( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._stressRatioThreshold = value + + @property + def criticalPorePressure( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the critical pore pressure value.""" + return self._criticalPorePressure + + @criticalPorePressure.setter + def criticalPorePressure( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._criticalPorePressure = value + + @property + def criticalPorePressureIndex( self: Self ) -> npt.NDArray[ np.float64 ] | None: + """Get the critical pore pressure index value.""" + return self._criticalPorePressureIndex + + @criticalPorePressureIndex.setter + def criticalPorePressureIndex( self: Self, value: npt.NDArray[ np.float64 ] ) -> None: + self._criticalPorePressureIndex = value + + def getAdvancedOutputValue( self: Self, name: str ) -> npt.NDArray[ np.float64 ] | None: + """Get the value of the advanced output wanted. + + Args: + name (str): The name of the advanced output. + + Returns: + npt.NDArray[ np.float64 ]: the value of the advanced output. + """ + if name == CRITICAL_TOTAL_STRESS_RATIO.attributeName: + return self.criticalTotalStressRatio + elif name == TOTAL_STRESS_RATIO_THRESHOLD.attributeName: + return self.stressRatioThreshold + elif name == CRITICAL_PORE_PRESSURE.attributeName: + return self.criticalPorePressure + elif name == CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName: + return self.criticalPorePressureIndex + else: + raise NameError + + physicalConstants: PhysicalConstants + _elasticModuli: ElasticModuliValue + _mandatoryAttributes: MandatoryAttributesValue + _basicOutput: BasicOutputValue + _advancedOutput: AdvancedOutputValue + + def __init__( + self: Self, + mesh: Union[ vtkPointSet, vtkUnstructuredGrid ], + computeAdvancedOutputs: bool = False, + speHandler: bool = False, + ) -> None: + """VTK Filter to perform Geomechanical output computation. + + Args: + mesh (Union[vtkPointsSet, vtkUnstructuredGrid]): Input mesh. + computeAdvancedOutputs (bool, optional): True to compute advanced geomechanical parameters, False otherwise. + Defaults to False. + speHandler (bool, optional): True to use a specific handler, False to use the internal handler. + Defaults to False. + """ + self.output: Union[ vtkPointSet, vtkUnstructuredGrid ] = mesh.NewInstance() + self.output.DeepCopy( mesh ) + + self.doComputeAdvancedOutputs: bool = computeAdvancedOutputs + self.physicalConstants = self.PhysicalConstants() + self._elasticModuli = self.ElasticModuliValue() + self._mandatoryAttributes = self.MandatoryAttributesValue() + self._basicOutput = self.BasicOutputValue() + self._advancedOutput = self.AdvancedOutputValue() + + self._attributesToCreate: list[ AttributeEnum ] = [] + + # Logger. + self.logger: Logger + if not speHandler: + self.logger = getLogger( loggerTitle, True ) + else: + self.logger = logging.getLogger( loggerTitle ) + self.logger.setLevel( logging.INFO ) + + def applyFilter( self: Self ) -> bool: + """Compute the geomechanical outputs of the mesh. + + Returns: + Boolean (bool): True if calculation successfully ended, False otherwise. + """ + if not self._checkMandatoryAttributes(): + return False + + if not self.computeBasicOutputs(): + return False + + if self.doComputeAdvancedOutputs and not self.computeAdvancedOutputs(): + return False + + # Create an attribute on the mesh for each geomechanics outputs computed: + for attribute in self._attributesToCreate: + attributeName: str = attribute.attributeName + onPoints: bool = attribute.isOnPoints + array: npt.NDArray[ np.float64 ] | None + if attribute in ELASTIC_MODULI: + array = self._elasticModuli.getElasticModulusValue( attributeName ) + elif attribute in BASIC_OUTPUTS: + array = self._basicOutput.getBasicOutputValue( attributeName ) + elif attribute in ADVANCED_OUTPUTS: + array = self._advancedOutput.getAdvancedOutputValue( attributeName ) + componentNames: tuple[ str, ...] = () + if attribute.nbComponent == 6: + componentNames = ComponentNameEnum.XYZ.value + + if not createAttribute( self.output, + array, + attributeName, + componentNames=componentNames, + onPoints=onPoints, + logger=self.logger ): + return False + + self.logger.info( "All the geomechanics properties have been added to the mesh." ) + self.logger.info( "The filter succeeded." ) + return True + + def getOutput( self: Self ) -> Union[ vtkPointSet, vtkUnstructuredGrid ]: + """Get the mesh with the geomechanical outputs as attributes. + + Returns: + Union[vtkPointSet, vtkUnstructuredGrid]: The mesh with the geomechanical outputs as attributes. + """ + return self.output + + def setLoggerHandler( self: Self, handler: logging.Handler ) -> None: + """Set a specific handler for the filter logger. + + In this filter 4 log levels are use, .info, .error, .warning and .critical, be sure to have at least the same 4 levels. + + Args: + handler (logging.Handler): The handler to add. + """ + if not self.logger.hasHandlers(): + self.logger.addHandler( handler ) + else: + self.logger.warning( + "The logger already has an handler, to use yours set the argument 'speHandler' to True during the filter initialization." + ) + + def getOutputType( self: Self ) -> str: + """Get output object type. + + Returns: + str: Type of output object. + """ + return self.output.GetClassName() + + def _checkMandatoryAttributes( self: Self ) -> bool: + """Check that mandatory attributes are present in the mesh. + + The mesh must contains: + - The Young modulus and the Poisson's ratio named "youngModulus" and "poissonRatio" or bulk and shear moduli named "bulkModulus" and "shearModulus" + - The initial Young modulus and Poisson's ratio named "youngModulusInitial" and "poissonRatioInitial" or the initial bulk modulus named "bulkModulusInitial" + - The porosity named "porosity" + - The initial porosity named "porosityInitial" + - The pressure named "pressure" + - The delta of pressure named "deltaPressure" + - The density named "density" + - The effective stress named "stressEffective" + - The initial effective stress named "stressEffectiveInitial" + + Returns: + bool: True if all needed attributes are present, False otherwise + """ + for elasticModulus in ELASTIC_MODULI: + elasticModulusName: str = elasticModulus.attributeName + elasticModulusOnPoints: bool = elasticModulus.isOnPoints + if isAttributeInObject( self.output, elasticModulusName, elasticModulusOnPoints ): + self._elasticModuli.setElasticModulusValue( + elasticModulus.attributeName, + getArrayInObject( self.output, elasticModulusName, elasticModulusOnPoints ) ) + + # Check the presence of the elastic moduli at the current time. + self.computeYoungPoisson: bool + if self._elasticModuli.youngModulus is None and self._elasticModuli.poissonRatio is None: + if self._elasticModuli.bulkModulus is not None and self._elasticModuli.shearModulus is not None: + self._elasticModuli.youngModulus = fcts.youngModulus( self._elasticModuli.bulkModulus, + self._elasticModuli.shearModulus ) + self._attributesToCreate.append( YOUNG_MODULUS ) + self._elasticModuli.poissonRatio = fcts.poissonRatio( self._elasticModuli.bulkModulus, + self._elasticModuli.shearModulus ) + self._attributesToCreate.append( POISSON_RATIO ) + self.computeYoungPoisson = True + else: + self.logger.error( + f"{ BULK_MODULUS.attributeName } or { SHEAR_MODULUS.attributeName } are missing to compute geomechanical outputs." + ) + return False + elif self._elasticModuli.bulkModulus is None and self._elasticModuli.shearModulus is None: + if self._elasticModuli.youngModulus is not None and self._elasticModuli.poissonRatio is not None: + self._elasticModuli.bulkModulus = fcts.bulkModulus( self._elasticModuli.youngModulus, + self._elasticModuli.poissonRatio ) + self._attributesToCreate.append( BULK_MODULUS ) + self._elasticModuli.shearModulus = fcts.shearModulus( self._elasticModuli.youngModulus, + self._elasticModuli.poissonRatio ) + self._attributesToCreate.append( SHEAR_MODULUS ) + self.computeYoungPoisson = False + else: + self.logger.error( + f"{ YOUNG_MODULUS.attributeName } or { POISSON_RATIO.attributeName } are missing to compute geomechanical outputs." + ) + return False + else: + self.logger.error( + f"{ BULK_MODULUS.attributeName } and { SHEAR_MODULUS.attributeName } or { YOUNG_MODULUS.attributeName } and { POISSON_RATIO.attributeName } are mandatory to compute geomechanical outputs." + ) + return False + + # Check the presence of the elastic moduli at the initial time. + if self._elasticModuli.bulkModulusT0 is None: + if self._elasticModuli.youngModulusT0 is not None and self._elasticModuli.poissonRatioT0 is not None: + self._elasticModuli.bulkModulusT0 = fcts.bulkModulus( self._elasticModuli.youngModulusT0, + self._elasticModuli.poissonRatioT0 ) + self._attributesToCreate.append( BULK_MODULUS_T0 ) + else: + self.logger.error( + f"{ BULK_MODULUS_T0.attributeName } or { YOUNG_MODULUS_T0.attributeName } and { POISSON_RATIO_T0.attributeName } are mandatory to compute geomechanical outputs." + ) + return False + + # Check the presence of the other mandatory attributes + for mandatoryAttribute in MANDATORY_ATTRIBUTES: + mandatoryAttributeName: str = mandatoryAttribute.attributeName + mandatoryAttributeOnPoints: bool = mandatoryAttribute.isOnPoints + if not isAttributeInObject( self.output, mandatoryAttributeName, mandatoryAttributeOnPoints ): + self.logger.error( + f"The mandatory attribute { mandatoryAttributeName } is missing to compute geomechanical outputs." ) + return False + else: + self._mandatoryAttributes.setMandatoryAttributeValue( + mandatoryAttributeName, + getArrayInObject( self.output, mandatoryAttributeName, mandatoryAttributeOnPoints ) ) + + return True + + def computeBasicOutputs( self: Self ) -> bool: + """Compute basic geomechanical outputs. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + if not self._computeBiotCoefficient(): + self.logger.error( "Biot coefficient computation failed." ) + return False + + if not self._computeCompressibilityCoefficient(): + return False + + if not self._computeRealEffectiveStressRatio(): + self.logger.error( "Effective stress ratio computation failed." ) + return False + + if not self._computeSpecificGravity(): + self.logger.error( "Specific gravity computation failed." ) + return False + + # TODO: deactivate lithostatic stress calculation until right formula + # if not self.computeLithostaticStress(): + # self.logger.error( "Lithostatic stress computation failed." ) + # return False + + if not self._computeTotalStresses(): + return False + + if not self._computeElasticStrain(): + self.logger.error( "Elastic strain computation failed." ) + return False + + if not self._computeEffectiveStressRatioOed(): + self.logger.error( "Effective stress ration in oedometric condition computation failed." ) + return False + + if not self._computeReservoirStressPathOed(): + self.logger.error( "Reservoir stress path in oedometric condition computation failed." ) + return False + + if not self._computeReservoirStressPathReal(): + return False + + self.logger.info( "All geomechanical basic outputs were successfully computed." ) + return True + + def computeAdvancedOutputs( self: Self ) -> bool: + """Compute advanced geomechanical outputs. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + if not self._computeCriticalTotalStressRatio(): + return False + + if not self._computeCriticalPorePressure(): + return False + + self.logger.info( "All geomechanical advanced outputs were successfully computed." ) + return True + + def _computeBiotCoefficient( self: Self ) -> bool: + """Compute Biot coefficient from default and grain bulk modulus. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + if not isAttributeInObject( self.output, BIOT_COEFFICIENT.attributeName, BIOT_COEFFICIENT.isOnPoints ): + self._basicOutput.biotCoefficient = fcts.biotCoefficient( self.physicalConstants.grainBulkModulus, + self._elasticModuli.bulkModulus ) + self._attributesToCreate.append( BIOT_COEFFICIENT ) + else: + self._basicOutput.biotCoefficient = getArrayInObject( self.output, BIOT_COEFFICIENT.attributeName, + BIOT_COEFFICIENT.isOnPoints ) + self.logger.warning( + f"{ BIOT_COEFFICIENT.attributeName } is already on the mesh, it has not been computed by the filter." ) + + return True + + def _computeCompressibilityCoefficient( self: Self ) -> bool: + """Compute compressibility coefficient from simulation outputs. + + Compressibility coefficient is computed from Poisson's ratio, bulk + modulus, Biot coefficient and Porosity. + + Returns: + bool: True if the attribute is correctly created, False otherwise. + """ + if not isAttributeInObject( self.output, COMPRESSIBILITY.attributeName, COMPRESSIBILITY.isOnPoints ): + self._basicOutput.compressibility = fcts.compressibility( self._elasticModuli.poissonRatio, + self._elasticModuli.bulkModulus, + self._basicOutput.biotCoefficient, + self._mandatoryAttributes.porosity ) + self._attributesToCreate.append( COMPRESSIBILITY ) + else: + self._basicOutput.compressibility = getArrayInObject( self.output, COMPRESSIBILITY.attributeName, + COMPRESSIBILITY.isOnPoints ) + self.logger.warning( + f"{ COMPRESSIBILITY.attributeName } is already on the mesh, it has not been computed by the filter." ) + + # oedometric compressibility + if not isAttributeInObject( self.output, COMPRESSIBILITY_OED.attributeName, COMPRESSIBILITY_OED.isOnPoints ): + self._basicOutput.compressibilityOed = fcts.compressibilityOed( self._elasticModuli.shearModulus, + self._elasticModuli.bulkModulus, + self._mandatoryAttributes.porosity ) + self._attributesToCreate.append( COMPRESSIBILITY_OED ) + else: + self._basicOutput.compressibilityOed = getArrayInObject( self.output, COMPRESSIBILITY_OED.attributeName, + COMPRESSIBILITY_OED.isOnPoints ) + self.logger.warning( + f"{ COMPRESSIBILITY_OED.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + # real compressibility + if not isAttributeInObject( self.output, COMPRESSIBILITY_REAL.attributeName, COMPRESSIBILITY_REAL.isOnPoints ): + self._basicOutput.compressibilityReal = fcts.compressibilityReal( + self._mandatoryAttributes.deltaPressure, self._mandatoryAttributes.porosity, + self._mandatoryAttributes.porosityInitial ) + self._attributesToCreate.append( COMPRESSIBILITY_REAL ) + else: + self._basicOutput.compressibilityReal = getArrayInObject( self.output, COMPRESSIBILITY_REAL.attributeName, + COMPRESSIBILITY_REAL.isOnPoints ) + self.logger.warning( + f"{ COMPRESSIBILITY_REAL.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + return True + + def _computeSpecificGravity( self: Self ) -> bool: + """Create Specific gravity attribute. + + Specific gravity is computed from rock density attribute and specific + density input. + + Returns: + bool: True if the attribute is correctly created, False otherwise. + """ + if not isAttributeInObject( self.output, SPECIFIC_GRAVITY.attributeName, SPECIFIC_GRAVITY.isOnPoints ): + self._basicOutput.specificGravity = fcts.specificGravity( self._mandatoryAttributes.density, + self.physicalConstants.specificDensity ) + self._attributesToCreate.append( SPECIFIC_GRAVITY ) + else: + self._basicOutput.specificGravity = getArrayInObject( self.output, SPECIFIC_GRAVITY.attributeName, + SPECIFIC_GRAVITY.isOnPoints ) + self.logger.warning( + f"{ SPECIFIC_GRAVITY.attributeName } is already on the mesh, it has not been computed by the filter." ) + + return True + + def _doComputeStressRatioReal( self: Self, stress: npt.NDArray[ np.float64 ], + basicOutput: AttributeEnum ) -> npt.NDArray[ np.float64 ]: + """Compute the ratio between horizontal and vertical effective stress. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + verticalStress: npt.NDArray[ np.float64 ] = stress[ :, 2 ] + # keep the minimum of the 2 horizontal components + horizontalStress: npt.NDArray[ np.float64 ] = np.min( stress[ :, :2 ], axis=1 ) + + stressRatioReal: npt.NDArray[ np.float64 ] + if not isAttributeInObject( self.output, basicOutput.attributeName, basicOutput.isOnPoints ): + stressRatioReal = fcts.stressRatio( horizontalStress, verticalStress ) + self._attributesToCreate.append( basicOutput ) + else: + stressRatioReal = getArrayInObject( self.output, basicOutput.attributeName, basicOutput.isOnPoints ) + self.logger.warning( + f"{ basicOutput.attributeName } is already on the mesh, it has not been computed by the filter." ) + + return stressRatioReal + + def _computeRealEffectiveStressRatio( self: Self ) -> bool: + """Compute effective stress ratio. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + if self._mandatoryAttributes.effectiveStress is not None: + self._basicOutput.effectiveStressRatioReal = self._doComputeStressRatioReal( + self._mandatoryAttributes.effectiveStress, STRESS_EFFECTIVE_RATIO_REAL ) + return True + else: + self.logger.error( + f"{ STRESS_EFFECTIVE_RATIO_REAL.attributeName } has not been computed, mandatory attribute { STRESS_EFFECTIVE.attributeName } is missing." + ) + return False + + def _doComputeTotalStress( + self: Self, + effectiveStress: npt.NDArray[ np.float64 ], + pressure: Union[ npt.NDArray[ np.float64 ], None ], + biotCoefficient: npt.NDArray[ np.float64 ], + ) -> npt.NDArray[ np.float64 ]: + """Compute total stress from effective stress and pressure. + + Args: + effectiveStress (npt.NDArray[np.float64]): Effective stress. + pressure (npt.NDArray[np.float64] | None): Pressure. + biotCoefficient (npt.NDArray[np.float64]): Biot coefficient. + + Returns: + npt.NDArray[ np.float64 ]: TotalStress. + """ + totalStress: npt.NDArray[ np.float64 ] + if pressure is None: + totalStress = np.copy( effectiveStress ) + self.logger.warning( "Pressure attribute is undefined, total stress will be equal to effective stress." ) + else: + if np.isnan( pressure ).any(): + self.logger.warning( "Some cells do not have pressure data, for those cells, pressure is set to 0." ) + pressure[ np.isnan( pressure ) ] = 0.0 + + totalStress = fcts.totalStress( effectiveStress, biotCoefficient, pressure ) + + return totalStress + + def _computeTotalStressInitial( self: Self ) -> bool: + """Compute total stress at initial time step. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + # Compute Biot at initial time step. + biotCoefficientT0: npt.NDArray[ np.float64 ] = fcts.biotCoefficient( self.physicalConstants.grainBulkModulus, + self._elasticModuli.bulkModulusT0 ) + + pressureT0: npt.NDArray[ np.float64 ] | None = None + # Case pressureT0 is None, total stress = effective stress + # (managed by doComputeTotalStress function) + if self._mandatoryAttributes.pressure is not None: + if self._mandatoryAttributes.deltaPressure is not None: + pressureT0 = self._mandatoryAttributes.pressure - self._mandatoryAttributes.deltaPressure + else: + self.logger.error( f"Mandatory attribute { DELTA_PRESSURE.attributeName } is missing." ) + return False + + if not isAttributeInObject( self.output, STRESS_TOTAL_T0.attributeName, STRESS_TOTAL_T0.isOnPoints ): + if self._mandatoryAttributes.effectiveStressT0 is not None: + self._basicOutput.totalStressT0 = self._doComputeTotalStress( + self._mandatoryAttributes.effectiveStressT0, pressureT0, biotCoefficientT0 ) + self._attributesToCreate.append( STRESS_TOTAL_T0 ) + else: + self.logger.error( + f"{ STRESS_TOTAL_T0.attributeName } has not been computed, mandatory attribute { STRESS_EFFECTIVE_T0.attributeName } is missing." + ) + return False + else: + self._basicOutput.totalStressT0 = getArrayInObject( self.output, STRESS_TOTAL_T0.attributeName, + STRESS_TOTAL_T0.isOnPoints ) + self.logger.warning( + f"{ STRESS_TOTAL_T0.attributeName } is already on the mesh, it has not been computed by the filter." ) + + return True + + def _computeTotalStresses( self: Self ) -> bool: + """Compute total stress total stress ratio. + + Total stress is computed at the initial and current time steps. + Total stress ratio is computed at current time step only. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + # Compute total stress at initial time step. + if not self._computeTotalStressInitial(): + self.logger.error( "Total stress at initial time step computation failed." ) + return False + + # Compute total stress at current time step. + if not isAttributeInObject( self.output, STRESS_TOTAL.attributeName, STRESS_TOTAL.isOnPoints ): + if self._mandatoryAttributes.effectiveStress is not None and self._basicOutput.biotCoefficient is not None: + self._basicOutput.totalStress = self._doComputeTotalStress( self._mandatoryAttributes.effectiveStress, + self._mandatoryAttributes.pressure, + self._basicOutput.biotCoefficient ) + self._attributesToCreate.append( STRESS_TOTAL ) + else: + self.logger.error( + f"{ STRESS_TOTAL.attributeName } has not been computed, mandatory attributes { STRESS_EFFECTIVE.attributeName } or { BIOT_COEFFICIENT.attributeName } are missing." + ) + return False + else: + self._basicOutput.totalStress = getArrayInObject( self.output, STRESS_TOTAL.attributeName, + STRESS_TOTAL.isOnPoints ) + self.logger.warning( + f"{ STRESS_TOTAL.attributeName } is already on the mesh, it has not been computed by the filter." ) + + # Compute total stress ratio. + if self._basicOutput.totalStress is not None: + self._basicOutput.totalStressRatioReal = self._doComputeStressRatioReal( self._basicOutput.totalStress, + STRESS_TOTAL_RATIO_REAL ) + else: + self.logger.error( + f"{ STRESS_TOTAL_RATIO_REAL.attributeName } has not been computed, Mandatory attribute { BIOT_COEFFICIENT.attributeName } is missing." + ) + return False + + return True + + def computeLithostaticStress( self: Self ) -> bool: + """Compute lithostatic stress. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + if not isAttributeInObject( self.output, LITHOSTATIC_STRESS.attributeName, LITHOSTATIC_STRESS.isOnPoints ): + depth: npt.NDArray[ np.float64 ] = self._doComputeDepthAlongLine( + ) if LITHOSTATIC_STRESS.isOnPoints else self._doComputeDepthInMesh() + self._basicOutput.lithostaticStress = fcts.lithostaticStress( depth, self._mandatoryAttributes.density, + GRAVITY ) + self._attributesToCreate.append( LITHOSTATIC_STRESS ) + else: + self._basicOutput.lithostaticStress = getArrayInObject( self.output, LITHOSTATIC_STRESS.attributeName, + LITHOSTATIC_STRESS.isOnPoints ) + self.logger.warning( + f"{ LITHOSTATIC_STRESS.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + return True + + def _doComputeDepthAlongLine( self: Self ) -> npt.NDArray[ np.float64 ]: + """Compute depth along a line. + + Returns: + npt.NDArray[np.float64]: 1D array with depth property + """ + # get z coordinate + zCoord: npt.NDArray[ np.float64 ] = self._getZcoordinates( True ) + assert zCoord is not None, "Depth coordinates cannot be computed." + + # TODO: to find how to compute depth in a general case + # GEOS z axis is upward + depth: npt.NDArray[ np.float64 ] = -1.0 * zCoord + return depth + + def _doComputeDepthInMesh( self: Self ) -> npt.NDArray[ np.float64 ]: + """Compute depth of each cell in a mesh. + + Returns: + npt.NDArray[np.float64]: array with depth property + """ + # get z coordinate + zCoord: npt.NDArray[ np.float64 ] = self._getZcoordinates( False ) + assert zCoord is not None, "Depth coordinates cannot be computed." + + # TODO: to find how to compute depth in a general case + depth: npt.NDArray[ np.float64 ] = -1.0 * zCoord + return depth + + def _getZcoordinates( self: Self, onPoints: bool ) -> npt.NDArray[ np.float64 ]: + """Get z coordinates from self.output. + + Args: + onPoints (bool): True if the attribute is on points, False if it is on cells. + + Returns: + npt.NDArray[np.float64]: 1D array with depth property + """ + # get z coordinate + zCoord: npt.NDArray[ np.float64 ] + pointCoords: npt.NDArray[ np.float64 ] = self._getPointCoordinates( onPoints ) + assert pointCoords is not None, "Point coordinates are undefined." + assert pointCoords.shape[ 1 ] == 2, "Point coordinates are undefined." + zCoord = pointCoords[ :, 2 ] + return zCoord + + def _getPointCoordinates( self: Self, onPoints: bool ) -> npt.NDArray[ np.float64 ]: + """Get the coordinates of Points or Cell center. + + Args: + onPoints (bool): True if the attribute is on points, False if it is on cells. + + Returns: + npt.NDArray[np.float64]: points/cell center coordinates + """ + if onPoints: + return self.output.GetPoints() # type: ignore[no-any-return] + else: + # Find cell centers + filter = vtkCellCenters() + filter.SetInputDataObject( self.output ) + filter.Update() + return filter.GetOutput().GetPoints() # type: ignore[no-any-return] + + def _computeElasticStrain( self: Self ) -> bool: + """Compute elastic strain from effective stress and elastic modulus. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + if self._mandatoryAttributes.effectiveStress is not None and self._mandatoryAttributes.effectiveStressT0 is not None: + deltaEffectiveStress = self._mandatoryAttributes.effectiveStress - self._mandatoryAttributes.effectiveStressT0 + + if not isAttributeInObject( self.output, STRAIN_ELASTIC.attributeName, STRAIN_ELASTIC.isOnPoints ): + if self.computeYoungPoisson: + self._basicOutput.elasticStrain = fcts.elasticStrainFromBulkShear( + deltaEffectiveStress, self._elasticModuli.bulkModulus, self._elasticModuli.shearModulus ) + else: + self._basicOutput.elasticStrain = fcts.elasticStrainFromYoungPoisson( + deltaEffectiveStress, self._elasticModuli.youngModulus, self._elasticModuli.poissonRatio ) + self._attributesToCreate.append( STRAIN_ELASTIC ) + else: + self._basicOutput.totalStressT0 = getArrayInObject( self.output, STRAIN_ELASTIC.attributeName, + STRAIN_ELASTIC.isOnPoints ) + self.logger.warning( + f"{ STRAIN_ELASTIC.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + return True + + else: + self.logger.error( + f"{ STRAIN_ELASTIC.attributeName } has not been computed, mandatory attributes { STRESS_EFFECTIVE.attributeName } or { STRESS_EFFECTIVE_T0.attributeName } are missing." + ) + return False + + def _computeReservoirStressPathReal( self: Self ) -> bool: + """Compute reservoir stress paths. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + # create delta stress attribute for QC + if not isAttributeInObject( self.output, STRESS_TOTAL_DELTA.attributeName, STRESS_TOTAL_DELTA.isOnPoints ): + if self._basicOutput.totalStress is not None and self._basicOutput.totalStressT0 is not None: + self._basicOutput.deltaTotalStress = self._basicOutput.totalStress - self._basicOutput.totalStressT0 + self._attributesToCreate.append( STRESS_TOTAL_DELTA ) + else: + self.logger.error( + f"{ STRESS_TOTAL_DELTA.attributeName } has not been computed, mandatory attributes { STRESS_TOTAL.attributeName } or { STRESS_TOTAL_T0.attributeName } are missing." + ) + return False + else: + self._basicOutput.deltaTotalStress = getArrayInObject( self.output, STRESS_TOTAL_DELTA.attributeName, + STRESS_TOTAL_DELTA.isOnPoints ) + self.logger.warning( + f"{ STRESS_TOTAL_DELTA.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + if not isAttributeInObject( self.output, RSP_REAL.attributeName, RSP_REAL.isOnPoints ): + self._basicOutput.rspReal = fcts.reservoirStressPathReal( self._basicOutput.deltaTotalStress, + self._mandatoryAttributes.deltaPressure ) + self._attributesToCreate.append( RSP_REAL ) + else: + self._basicOutput.rspReal = getArrayInObject( self.output, RSP_REAL.attributeName, RSP_REAL.isOnPoints ) + self.logger.warning( + f"{ RSP_REAL.attributeName } is already on the mesh, it has not been computed by the filter." ) + + return True + + def _computeReservoirStressPathOed( self: Self ) -> bool: + """Compute Reservoir Stress Path in oedometric conditions. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + if not isAttributeInObject( self.output, RSP_OED.attributeName, RSP_OED.isOnPoints ): + self._basicOutput.rspOed = fcts.reservoirStressPathOed( self._basicOutput.biotCoefficient, + self._elasticModuli.poissonRatio ) + self._attributesToCreate.append( RSP_OED ) + else: + self._basicOutput.rspOed = getArrayInObject( self.output, RSP_OED.attributeName, RSP_OED.isOnPoints ) + self.logger.warning( + f"{ RSP_OED.attributeName } is already on the mesh, it has not been computed by the filter." ) + + return True + + def _computeEffectiveStressRatioOed( self: Self ) -> bool: + """Compute the effective stress ratio in oedometric conditions. + + Returns: + bool: True if calculation successfully ended, False otherwise. + """ + if not isAttributeInObject( self.output, STRESS_EFFECTIVE_RATIO_OED.attributeName, + STRESS_EFFECTIVE_RATIO_OED.isOnPoints ): + self._basicOutput.effectiveStressRatioOed = fcts.deviatoricStressPathOed( self._elasticModuli.poissonRatio ) + self._attributesToCreate.append( STRESS_EFFECTIVE_RATIO_OED ) + else: + self._basicOutput.effectiveStressRatioOed = getArrayInObject( self.output, + STRESS_EFFECTIVE_RATIO_OED.attributeName, + STRESS_EFFECTIVE_RATIO_OED.isOnPoints ) + self.logger.warning( + f"{ STRESS_EFFECTIVE_RATIO_OED.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + return True + + def _computeCriticalTotalStressRatio( self: Self ) -> bool: + """Compute fracture index and fracture threshold. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + if not isAttributeInObject( self.output, CRITICAL_TOTAL_STRESS_RATIO.attributeName, + CRITICAL_TOTAL_STRESS_RATIO.isOnPoints ): + if self._basicOutput.totalStress is not None: + verticalStress: npt.NDArray[ np.float64 ] = self._basicOutput.totalStress[ :, 2 ] + self._advancedOutput.criticalTotalStressRatio = fcts.criticalTotalStressRatio( + self._mandatoryAttributes.pressure, verticalStress ) + self._attributesToCreate.append( CRITICAL_TOTAL_STRESS_RATIO ) + else: + self.logger.error( + f"{ CRITICAL_TOTAL_STRESS_RATIO.attributeName } has not been computed, mandatory attribute { STRESS_TOTAL.attributeName } is missing." + ) + return False + else: + self._advancedOutput.criticalTotalStressRatio = getArrayInObject( self.output, + CRITICAL_TOTAL_STRESS_RATIO.attributeName, + CRITICAL_TOTAL_STRESS_RATIO.isOnPoints ) + self.logger.warning( + f"{ CRITICAL_TOTAL_STRESS_RATIO.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + if not isAttributeInObject( self.output, TOTAL_STRESS_RATIO_THRESHOLD.attributeName, + TOTAL_STRESS_RATIO_THRESHOLD.isOnPoints ): + if self._basicOutput.totalStress is not None: + mask: npt.NDArray[ np.bool_ ] = np.argmin( np.abs( self._basicOutput.totalStress[ :, :2 ] ), axis=1 ) + horizontalStress: npt.NDArray[ np.float64 ] = self._basicOutput.totalStress[ :, :2 ][ + np.arange( self._basicOutput.totalStress[ :, :2 ].shape[ 0 ] ), mask ] + self._advancedOutput.stressRatioThreshold = fcts.totalStressRatioThreshold( + self._mandatoryAttributes.pressure, horizontalStress ) + self._attributesToCreate.append( TOTAL_STRESS_RATIO_THRESHOLD ) + else: + self.logger.error( + f"{ TOTAL_STRESS_RATIO_THRESHOLD.attributeName } has not been computed, mandatory attribute { STRESS_TOTAL.attributeName } is missing." + ) + return False + else: + self._advancedOutput.stressRatioThreshold = getArrayInObject( self.output, + TOTAL_STRESS_RATIO_THRESHOLD.attributeName, + TOTAL_STRESS_RATIO_THRESHOLD.isOnPoints ) + self.logger.warning( + f"{ TOTAL_STRESS_RATIO_THRESHOLD.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + return True + + def _computeCriticalPorePressure( self: Self ) -> bool: + """Compute the critical pore pressure and the pressure index. + + Returns: + bool: return True if calculation successfully ended, False otherwise. + """ + if not isAttributeInObject( self.output, CRITICAL_PORE_PRESSURE.attributeName, + CRITICAL_PORE_PRESSURE.isOnPoints ): + if self._basicOutput.totalStress is not None: + self._advancedOutput.criticalPorePressure = fcts.criticalPorePressure( + -1.0 * self._basicOutput.totalStress, self.physicalConstants.rockCohesion, + self.physicalConstants.frictionAngle ) + self._attributesToCreate.append( CRITICAL_PORE_PRESSURE ) + else: + self.logger.error( + f"{ CRITICAL_PORE_PRESSURE.attributeName } has not been computed, mandatory attribute { STRESS_TOTAL.attributeName } is missing." + ) + return False + else: + self._advancedOutput.criticalPorePressure = getArrayInObject( self.output, + CRITICAL_PORE_PRESSURE.attributeName, + CRITICAL_PORE_PRESSURE.isOnPoints ) + self.logger.warning( + f"{ CRITICAL_PORE_PRESSURE.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + # Add critical pore pressure index (i.e., ratio between pressure and criticalPorePressure) + if not isAttributeInObject( self.output, CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName, + CRITICAL_PORE_PRESSURE_THRESHOLD.isOnPoints ): + self._advancedOutput.criticalPorePressureIndex = fcts.criticalPorePressureThreshold( + self._mandatoryAttributes.pressure, self._advancedOutput.criticalPorePressure ) + self._attributesToCreate.append( CRITICAL_PORE_PRESSURE_THRESHOLD ) + else: + self._advancedOutput.criticalPorePressureIndex = getArrayInObject( + self.output, CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName, + CRITICAL_PORE_PRESSURE_THRESHOLD.isOnPoints ) + self.logger.warning( + f"{ CRITICAL_PORE_PRESSURE_THRESHOLD.attributeName } is already on the mesh, it has not been computed by the filter." + ) + + return True diff --git a/geos-pv/requirements.txt b/geos-pv/requirements.txt index edb1046c..6ec1b5a9 100644 --- a/geos-pv/requirements.txt +++ b/geos-pv/requirements.txt @@ -1,4 +1,5 @@ geos-geomechanics geos-mesh geos-posp -geos-utils \ No newline at end of file +geos-utils +geos-processing \ No newline at end of file diff --git a/geos-pv/src/geos/pv/plugins/PVGeomechanicsCalculator.py b/geos-pv/src/geos/pv/plugins/PVGeomechanicsCalculator.py new file mode 100644 index 00000000..b6f0afff --- /dev/null +++ b/geos-pv/src/geos/pv/plugins/PVGeomechanicsCalculator.py @@ -0,0 +1,286 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Martin Lemay, Romain Baville +# ruff: noqa: E402 # disable Module level import not at top of file +import sys +from pathlib import Path +from typing import Union +from typing_extensions import Self + +from paraview.util.vtkAlgorithm import ( # type: ignore[import-not-found] + VTKPythonAlgorithmBase, smdomain, smhint, smproperty, smproxy, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/util/vtkAlgorithm.py +from paraview.detail.loghandler import ( # type: ignore[import-not-found] + VTKHandler, +) # source: https://github.com/Kitware/ParaView/blob/master/Wrapping/Python/paraview/detail/loghandler.py + +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector +from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid + +# update sys.path to load all GEOS Python Package dependencies +geos_pv_path: Path = Path( __file__ ).parent.parent.parent.parent.parent +sys.path.insert( 0, str( geos_pv_path / "src" ) ) +from geos.pv.utils.config import update_paths + +update_paths() + +from geos.utils.PhysicalConstants import ( + DEFAULT_FRICTION_ANGLE_RAD, + DEFAULT_GRAIN_BULK_MODULUS, + DEFAULT_ROCK_COHESION, + WATER_DENSITY, +) +from geos.processing.post_processing.GeomechanicsCalculator import GeomechanicsCalculator + +__doc__ = """ +PVGeomechanicsCalculator is a paraview plugin that allows to compute additional +Geomechanical properties from existing ones. + +The basic geomechanics outputs are: + - Elastic modulus (young modulus and poisson ratio or bulk modulus and shear modulus) + - Biot coefficient + - Compressibility, oedometric compressibility and real compressibility coefficient + - Specific gravity + - Real effective stress ratio + - Total initial stress, total current stress and total stress ratio + - Lithostatic stress (physics to update) + - Elastic stain + - Reservoir stress path and reservoir stress path in oedometric condition + +The advanced geomechanics outputs are: + - fracture index and threshold + - Critical pore pressure and pressure index + +PVGeomechanicsCalculator paraview plugin input mesh is either vtkPointSet or vtkUnstructuredGrid +and returned mesh is of same type as input. + +.. Important:: + Please refer to the PVGeosExtractMergeBlockVolume* plugins to provide the correct input. + +To use it: + +* Load the module in Paraview: Tools > Manage Plugins... > Load new > PVGeomechanicsCalculator +* Select the mesh you want to compute geomechanics output on +* Search Filters > 3- Geos Geomechanics > Geos Geomechanics Calculator +* Set physical constants and computeAdvancedOutput if needed +* Apply + +""" + + +@smproxy.filter( name="PVGeomechanicsCalculator", label="Geos Geomechanics Calculator" ) +@smhint.xml( """""" ) +@smproperty.input( name="Input", port_index=0 ) +@smdomain.datatype( dataTypes=[ "vtkUnstructuredGrid", "vtkPointSet" ], composite_data_supported=True ) +class PVGeomechanicsCalculator( VTKPythonAlgorithmBase ): + + def __init__( self: Self ) -> None: + """Paraview plugin to compute additional geomechanical outputs. + + Input is either a vtkUnstructuredGrid or vtkPointSet with Geos + geomechanical properties. + """ + super().__init__( nInputPorts=1, nOutputPorts=1, outputType="vtkPointSet" ) + + self.computeAdvancedOutputs: bool = False + + # Defaults physical constants + ## Basic outputs + self.grainBulkModulus: float = DEFAULT_GRAIN_BULK_MODULUS + self.specificDensity: float = WATER_DENSITY + ## Advanced outputs + self.rockCohesion: float = DEFAULT_ROCK_COHESION + self.frictionAngle: float = DEFAULT_FRICTION_ANGLE_RAD + + @smproperty.doublevector( + name="GrainBulkModulus", + label="Grain bulk modulus (Pa)", + default_values=DEFAULT_GRAIN_BULK_MODULUS, + panel_visibility="default", + ) + @smdomain.xml( """ + + Reference grain bulk modulus to compute Biot coefficient. + The unit is Pa. Default is Quartz bulk modulus (i.e., 38GPa). + + """ ) + def setGrainBulkModulus( self: Self, grainBulkModulus: float ) -> None: + """Set grain bulk modulus. + + Args: + grainBulkModulus (float): Grain bulk modulus (Pa). + """ + self.grainBulkModulus = grainBulkModulus + self.Modified() + + @smproperty.doublevector( + name="SpecificDensity", + label="Specific Density (kg/m3)", + default_values=WATER_DENSITY, + panel_visibility="default", + ) + @smdomain.xml( """ + + Reference density to compute specific gravity. + The unit is kg/m3. Default is fresh water density (i.e., 1000 kg/m3). + + """ ) + def setSpecificDensity( self: Self, specificDensity: float ) -> None: + """Set specific density. + + Args: + specificDensity (float): Reference specific density (kg/m3). + """ + self.specificDensity = specificDensity + self.Modified() + + @smproperty.xml( """ + + + + + """ ) + def groupBasicOutputParameters( self: Self ) -> None: + """Organize groups.""" + self.Modified() + + @smproperty.intvector( + name="ComputeAdvancedOutputs", + label="Compute advanced geomechanical outputs", + default_values=0, + panel_visibility="default", + ) + @smdomain.xml( """ + + + Check to compute advanced geomechanical outputs including + reservoir stress paths and fracture indexes. + + """ ) + def setComputeAdvancedOutputs( self: Self, computeAdvancedOutputs: bool ) -> None: + """Set advanced output calculation option. + + Args: + computeAdvancedOutputs (bool): True to compute advanced geomechanical parameters, False otherwise. + """ + self.computeAdvancedOutputs = computeAdvancedOutputs + self.Modified() + + @smproperty.doublevector( + name="RockCohesion", + label="Rock Cohesion (Pa)", + default_values=DEFAULT_ROCK_COHESION, + panel_visibility="default", + ) + @smdomain.xml( """ + + Reference rock cohesion to compute critical pore pressure. + The unit is Pa.Default is fractured case (i.e., 0. Pa). + + """ ) + def setRockCohesion( self: Self, rockCohesion: float ) -> None: + """Set rock cohesion. + + Args: + rockCohesion (float): Rock cohesion (Pa). + """ + self.rockCohesion = rockCohesion + self.Modified() + + @smproperty.doublevector( + name="FrictionAngle", + label="Friction Angle (rad)", + default_values=DEFAULT_FRICTION_ANGLE_RAD, + panel_visibility="default", + ) + @smdomain.xml( """ + + Reference friction angle to compute critical pore pressure. + The unit is rad. Default is 10.0 / 180.0 * np.pi rad. + + """ ) + def setFrictionAngle( self: Self, frictionAngle: float ) -> None: + """Set friction angle. + + Args: + frictionAngle (float): Friction angle (rad). + """ + self.frictionAngle = frictionAngle + self.Modified() + + @smproperty.xml( """ + + + + + + + + """ ) + def groupAdvancedOutputParameters( self: Self ) -> None: + """Organize groups.""" + self.Modified() + + def RequestDataObject( + self: Self, + request: vtkInformation, + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestDataObject. + + Args: + request (vtkInformation): Request. + inInfoVec (list[vtkInformationVector]): Input objects. + outInfoVec (vtkInformationVector): Output objects. + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inData = self.GetInputData( inInfoVec, 0, 0 ) + outData = self.GetOutputData( outInfoVec, 0 ) + assert inData is not None + if outData is None or ( not outData.IsA( inData.GetClassName() ) ): + outData = inData.NewInstance() + outInfoVec.GetInformationObject( 0 ).Set( outData.DATA_OBJECT(), outData ) + return super().RequestDataObject( request, inInfoVec, outInfoVec ) # type: ignore[no-any-return] + + def RequestData( + self: Self, + request: vtkInformation, # noqa: F841 + inInfoVec: list[ vtkInformationVector ], + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): Request. + inInfoVec (list[vtkInformationVector]): Input objects. + outInfoVec (vtkInformationVector): Output objects. + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + inputMesh: Union[ vtkPointSet, vtkUnstructuredGrid ] = self.GetInputData( inInfoVec, 0, 0 ) + outputMesh: Union[ vtkPointSet, vtkUnstructuredGrid ] = self.GetOutputData( outInfoVec, 0 ) + assert inputMesh is not None, "Input server mesh is null." + assert outputMesh is not None, "Output pipeline is null." + + filter: GeomechanicsCalculator = GeomechanicsCalculator( inputMesh, self.computeAdvancedOutputs, True ) + + if not filter.logger.hasHandlers(): + filter.setLoggerHandler( VTKHandler() ) + + filter.physicalConstants.grainBulkModulus = self.grainBulkModulus + filter.physicalConstants.specificDensity = self.specificDensity + filter.physicalConstants.rockCohesion = self.rockCohesion + filter.physicalConstants.frictionAngle = self.frictionAngle + + if filter.applyFilter(): + outputMesh.ShallowCopy( filter.getOutput() ) + outputMesh.Modified() + + return 1 diff --git a/geos-utils/src/geos/utils/GeosOutputsConstants.py b/geos-utils/src/geos/utils/GeosOutputsConstants.py index 3e684aaa..8a53dca1 100644 --- a/geos-utils/src/geos/utils/GeosOutputsConstants.py +++ b/geos-utils/src/geos/utils/GeosOutputsConstants.py @@ -6,7 +6,7 @@ from typing_extensions import Self __doc__ = """ -GeosOutputsConstants module defines usefull constant names such as attribute +GeosOutputsConstants module defines useful constant names such as attribute names, domain names, phase types, and the lists of attribute names to process. .. WARNING:: @@ -28,7 +28,7 @@ class AttributeEnum( Enum ): def __init__( self: Self, attributeName: str, nbComponent: int, onPoints: bool ) -> None: - """Define the enumeration to store attrbute properties. + """Define the enumeration to store attribute properties. Args: attributeName (str): name of the attribute @@ -134,7 +134,7 @@ class GeosMeshOutputsEnum( AttributeEnum ): DELTA_PRESSURE = ( "deltaPressure", 1, False ) MASS = ( "mass", 1, False ) - # geomechanic attributes + # geomechanics attributes ROCK_DENSITY = ( "density", 1, False ) PERMEABILITY = ( "permeability", 1, False ) POROSITY = ( "porosity", 1, False ) @@ -182,11 +182,11 @@ class PostProcessingOutputsEnum( AttributeEnum ): SPECIFIC_GRAVITY = ( "specificGravity", 1, False ) LITHOSTATIC_STRESS = ( "stressLithostatic", 1, False ) STRESS_EFFECTIVE_INITIAL = ( "stressEffectiveInitial", 6, False ) - STRESS_EFFECTIVE_RATIO_REAL = ( "stressEffectiveRatio_real", 6, False ) - STRESS_EFFECTIVE_RATIO_OED = ( "stressEffectiveRatio_oed", 6, False ) + STRESS_EFFECTIVE_RATIO_REAL = ( "stressEffectiveRatio_real", 1, False ) + STRESS_EFFECTIVE_RATIO_OED = ( "stressEffectiveRatio_oed", 1, False ) STRESS_TOTAL = ( "stressTotal", 6, False ) STRESS_TOTAL_INITIAL = ( "stressTotalInitial", 6, False ) - STRESS_TOTAL_RATIO_REAL = ( "stressTotalRatio_real", 6, False ) + STRESS_TOTAL_RATIO_REAL = ( "stressTotalRatio_real", 1, False ) STRESS_TOTAL_DELTA = ( "deltaStressTotal", 6, False ) STRAIN_ELASTIC = ( "strainElastic", 6, False ) RSP_OED = ( "rsp_oed", 1, False ) diff --git a/install_packages.sh b/install_packages.sh index 74238daf..e7160228 100755 --- a/install_packages.sh +++ b/install_packages.sh @@ -3,6 +3,7 @@ python -m pip install --upgrade ./geos-utils python -m pip install --upgrade ./geos-geomechanics python -m pip install --upgrade ./geos-mesh python -m pip install --upgrade ./geos-posp +python -m pip install --upgrade ./geos-processing python -m pip install --upgrade ./geos-xml-tools python -m pip install --upgrade ./geos-xml-viewer python -m pip install --upgrade ./hdf5-wrapper