diff --git a/pygeos-tools/pyproject.toml b/pygeos-tools/pyproject.toml index 24f846fd..2a6f5dec 100644 --- a/pygeos-tools/pyproject.toml +++ b/pygeos-tools/pyproject.toml @@ -22,10 +22,6 @@ dependencies = [ "numpy", "scipy", "mpi4py", - "vtk", - "pyevtk", - "xmltodict", - "h5py", ] [tool.mypy] diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py b/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py deleted file mode 100644 index f6faf41b..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py +++ /dev/null @@ -1,202 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import sys -import argparse - -class GeosxArgs: - """ - Class containing the main argument in command line for GEOSX - - Attributes - ---------- - cmdline : list of str - List corresponding to a splitted GEOSX submission line - options : dict - Dict containing GEOSX main options - """ - def __init__(self, args=sys.argv.copy()): - """ - Parameters - ---------- - args : list of str - List corresponding to a splitted submission line with options - """ - self.cmdline = args - self.options = self.optionsParser() - - - def optionsParser(self, cmdline=[]): - """ - Return a dict with useful geosx options parsed from a list/submission line. - - Parameters - ---------- - cmdline : list of str - List corresponding to a splitted GEOSX submission line - - Returns - -------- - dict : - Dict containing GEOSX main options - """ - if not cmdline: - cmdline = self.cmdline - desc = "Parser of GEOSX command line" - parser = argparse.ArgumentParser("GEOSX command line parser", - allow_abbrev=False, - description=desc) - - parser.add_argument("--input", "-i", "--xml", - type=str, dest="xml", default=None, - help="XML input file") - parser.add_argument("--restart", "-r", - dest="restart", type=str, default=None, - help="Target restart filename") - parser.add_argument("-x-partitions", "-x", - type=int, dest="xpart", default=None, - help="Number of partitions in the x direction.") - parser.add_argument("-y-partitions", "-y", - type=int, dest="ypart", default=None, - help="Number of partitions in the y direction.") - parser.add_argument("-z-partitions", "-z", - type=int, dest="zpart", default=None, - help="Number of partitions in the z direction.") - parser.add_argument("--output", "-o", - type=str, dest="outputdir", default=None, - help="Directory to put the output files") - parser.add_argument("--use-nonblocking", "-b", default=None, - dest="non_blocking", action="store_true", - help="Use non-blocking MPI communication") - parser.add_argument("--name", "-n", - dest="name", default=None, - help="Name of the problem, used for output") - parser.add_argument("--suppress-pinned", "-s", - dest="supp_pinned", action="store_true", default=None, - help="Suppress usage of pinned memory for MPI communication buffers") - parser.add_argument("--timers", "-t", - type=str, dest="timers", default=None, - help="String specifying the type of timer output") - parser.add_argument("--trace-data-migration", - dest="trace_data_mig", action="store_true", default=None, - help="Trace host-device data migration") - parser.add_argument("--pause-for", - dest="pause_for", default=None, - help="Pause geosx for a given number of seconds before starting execution") - - return vars(parser.parse_known_args(cmdline)[0]) - - - def updateArg(self, optionName=None, newValue=None): - """ - Update the GEOSX initialization arguments - - Parameters - ---------- - optionName : str - Name of the option to update - newValue : str - New value for the argument to be updated. - - Returns - ------- - bool : True if the option has been (found and) updated. False otherwise. - """ - if optionName is not None: - if self.options[optionName] != newValue: - self.options.update({optionName: newValue}) - return True - - return False - - - def getCommandLine(self): - """ - Return the command line specific to GEOSX initialization - - Returns - -------- - cl : list of str - List containing all the options requested - """ - cl = [""] - for opt, val in self.options.items(): - if val is not None: - ab = GeosxAbbrevOption().getAbbrv(opt) - cl += [ab] - if not isinstance(self.options[opt], bool): - cl += [str(self.options[opt])] - - return cl - - -class GeosxAbbrevOption: - """ - Class containing GEOSX command line options abbreviations - - Attributes - ---------- - abbrvDict : dict - Dict containing lists of abbreviations for GEOSX command line options - """ - def __init__(self): - self.abbrvDict = {"xml": ("--input", "-i"), - "restart": ("--restart", "r"), - "xpart": ("-x-partitions", "-x"), - "ypart": ("-y-partitions", "-y"), - "zpart": ("-z-partitions", "-z"), - "non_blocking": ("--use-non-blocking", "-b"), - "name": ("--name", "-n"), - "supp_pinned": ("--suppress-pinned"), - "outputdir": ("--output", "-o"), - "timers": ("--timers", "-t"), - "trace_data_mig": ("--trace-data-migration"), - "pause_for": ("--pause-for") - } - - def getAbbrv(self, optionName=None): - """ - Returns the abbreviation corresponding to the requested option - - Parameters - ---------- - optionName : str - Name of the requested option - - Returns - ------- - str : Requested abbreviations - """ - return self.abbrvDict[optionName][-1] - - - def getAllAbbrv(self, optionName=None): - """ - Returns the abbreviation corresponding to the requested option - - Parameters - ---------- - optionName : str - Name of the requested option - - Returns - ------- - list of str : - Requested abbreviations - """ - try: - return self.abbrvDict[optionName] - - except: - return [""] diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py b/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py deleted file mode 100644 index 161fdab7..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py +++ /dev/null @@ -1,120 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import os -from xml.etree import cElementTree as ET -import xmltodict -from re import findall -from ..mesh.InternalMesh import InternalMesh -from ..mesh.VtkMesh import VTKMesh - - -class XML(): - def __init__(self, xmlFile): - self.filename = xmlFile - - self.tree = ET.parse(xmlFile) - root = self.tree.getroot() - to_string = ET.tostring(root, method='xml') - - self.outputs = None - - root = xmltodict.parse(to_string, attr_prefix="" ,dict_constructor=dict) - for k, v in root['Problem'].items(): - words = findall('[A-Z][^A-Z]*', k) - words[0] = words[0].lower() - attr = "" - for word in words: - attr += word - setattr(self, attr, v) - - - - def updateSolvers(self, solverName, **kwargs): - root = self.tree.getroot() - solver = root.find("./Solvers/"+solverName) - for k, v in kwargs.items(): - if k in solver.attrib: - solver.set(k, v) - self.solvers[solverName].update({k:str(v)}) - - - def updateMesh(self, **kwargs): - root = self.tree.getroot() - mesh = root.find("./Mesh//") - for k, v in kwargs.items(): - if k in mesh.attrib: - mesh.set(k, v) - self.mesh[mesh.tag].update({k:str(v)}) - - - def updateGeometry(self, boxname, **kwargs): - root = self.tree.getroot() - geometry = root.find("./Geometry//*[@name="+boxname+"]") - - for i in len(self.geometry[geometry.tag]): - box = self.geometry[geometry.tag][i] - if boxname == box["name"]: - break - - for k, v in kwargs.items(): - if k in geometry.attrib: - geometry.set(k, v) - self.geometry[geometry.tag][i].update({k:str(v)}) - - - def getMeshObject(self): - - if "InternalMesh" in self.mesh.keys(): - #Not working properly for now - return InternalMesh(self) - - elif "VTKMesh" in self.mesh.keys(): - vtkFile = self.mesh["VTKMesh"]["file"] - if not os.path.isabs(vtkFile): - vtkFile = os.path.join(os.path.split(self.filename)[0], vtkFile) - return VTKMesh(vtkFile) - - def getAttribute(self, parentElement, attributeTag): - if parentElement == "root": - pElement = self.tree.find(f"./[@{attributeTag}]") - else: - pElement = self.tree.find(f"./*/{parentElement}/[@{attributeTag}]") - - return pElement.get(attributeTag) - - - def getSolverType(self): - return [k for k in self.solvers.keys() if k[0].isupper()] - - - def getSourcesAndReceivers(self): - solverType = self.getSolverType() - if len(solverType) > 1: - pass - else: - src = self.getAttribute(f"{solverType[0]}", "sourceCoordinates") - src = eval(src.replace("{", "[").replace("}", "]")) - - rcv = self.getAttribute(f"{solverType[0]}", "receiverCoordinates") - rcv = eval(rcv.replace("{", "[").replace("}", "]")) - - return src, rcv - - - def exportToXml(self, filename=None): - if filename is None: - self.tree.write(self.filename) - else: - self.tree.write(filename) diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py b/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py deleted file mode 100644 index af3b85b4..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py +++ /dev/null @@ -1,18 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -"""Input handle""" - -from .GeosxArgs import GeosxArgs, GeosxAbbrevOption -from .Xml import XML \ No newline at end of file diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py deleted file mode 100644 index 5019b913..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py +++ /dev/null @@ -1,144 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import numpy as np - - -class InternalMesh: - """ - GEOSX Internal Mesh - - Attributes - ---------- - xml : XML - XML object containing the information on the mesh - bounds : list of list - Real bounds of the mesh [[xmin, xmax],[ymin,ymax],[zmin,zmax]] - nx : int - Number of elements in the x direction - ny : int - Number of elements in the y direction - nz : int - Number of elements in the z direction - order : int - Mesh order - cellBlockNames : str - Names of each mesh block - cellBounds : - elementTypes : - Element types of each mesh block - numberOfCells : int - Total number of cells - numberOfPoints : int - Total number of points - fieldSpecifications : dict - Dict containing the mesh field specifications - """ - def __init__(self, xml): - """ - Parameters - ---------- - xml : XML - XML object containing the information on the mesh - """ - self.xml = xml - - mesh = xml.mesh["InternalMesh"] - elementRegion = xml.elementRegions["CellElementRegion"] - fieldSpecifications = xml.fieldSpecifications - - name = mesh["name"] - - xCoords = list(eval(mesh["xCoords"])) - yCoords = list(eval(mesh["yCoords"])) - zCoords = list(eval(mesh["zCoords"])) - - self.bounds = [[xCoords[0], xCoords[-1]], [yCoords[0], yCoords[-1]], [zCoords[0], zCoords[-1]]] - - - nxStr = mesh["nx"].strip('').replace('{','').replace('}','').split(',') - nyStr = mesh["ny"].strip('').replace('{','').replace('}','').split(',') - nzStr = mesh["nz"].strip('').replace('{','').replace('}','').split(',') - - nx = [eval(nx) for nx in nxStr] - ny = [eval(ny) for ny in nyStr] - nz = [eval(nz) for nz in nzStr] - - self.nx = nx - self.ny = ny - self.nz = nz - - order = 1 - self.order = order - - self.cellBlockNames = mesh["cellBlockNames"].strip('').replace('{','').replace('}','').split(',') - - xlayers = [] - ylayers = [] - zlayers = [] - for i in range(len(nx)): - xlayers.append([xCoords[i], xCoords[i+1]]) - for i in range(len(ny)): - ylayers.append([yCoords[i], yCoords[i+1]]) - for i in range(len(nz)): - zlayers.append([zCoords[i], zCoords[i+1]]) - - self.layers = [xlayers, ylayers, zlayers] - - xCellsBounds = np.zeros(sum(nx)+1) - yCellsBounds = np.zeros(sum(ny)+1) - zCellsBounds = np.zeros(sum(nz)+1) - - for i in range(len(nx)): - xstep = (xlayers[i][1]-xlayers[i][0])/nx[i] - if i == 0: - xCellsBounds[0:nx[i]] = np.arange(xlayers[i][0], xlayers[i][1], xstep) - else : - xCellsBounds[nx[i-1]:sum(nx[0:i+1])] = np.arange(xlayers[i][0], xlayers[i][1], xstep) - xCellsBounds[nx[-1]] = xlayers[i][1] - - for i in range(len(ny)): - ystep = (ylayers[i][1]-ylayers[i][0])/ny[i] - if i == 0: - yCellsBounds[0:ny[i]] = np.arange(ylayers[i][0], ylayers[i][1], ystep) - else : - xCellsBounds[ny[i-1]:sum(ny[0:i+1])] = np.arange(ylayers[i][0], ylayers[i][1], ystep) - yCellsBounds[ny[-1]] = ylayers[i][1] - - for i in range(len(nz)): - zstep = (zlayers[i][1]-zlayers[i][0])/nz[i] - if i == 0: - zCellsBounds[0:nz[i]] = np.arange(zlayers[i][0], zlayers[i][1], zstep) - else : - zCellsBounds[nz[i-1]:sum(nz[0:i+1])] = np.arange(zlayers[i][0], zlayers[i][1], zstep) - zCellsBounds[nz[-1]] = zlayers[i][1] - - self.cellBounds = [xCellsBounds, yCellsBounds, zCellsBounds] - - elementTypes = mesh["elementTypes"].strip('').replace('{','').replace('}','').split(',') - - self.elementTypes=[] - for type in elementTypes: - if type == "C3D8": - self.elementTypes.append("Hexahedron") - else: - self.elementTypes.append(type) - - - self.numberOfCells = sum(nx) * sum(ny) * sum(nz) - self.numberOfPoints = (sum(nx) + 1) * (sum(ny) + 1) * (sum(nz) + 1) - - self.fieldSpecifications = xml.fieldSpecifications - self.isSet = True - diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py deleted file mode 100644 index cba7259e..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py +++ /dev/null @@ -1,197 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import copy -import numpy as np -from vtk.util import numpy_support as VN - - -class VTKFieldSpecifications: - """ - Object containing field dataset of a VTK format mesh - - Attributes - ---------- - arrays : dict - Dict containing the {array names : array} of the given dataset - """ - def __init__(self, fieldData): - """ - Parameters - ---------- - fieldData : vtk.vtkFieldData - Data contained in the VTK mesh - """ - self.arrays = {} - assert(fieldData.IsA("vtkFieldData")) - for i in range(fieldData.GetNumberOfArrays()): - value = fieldData.GetArray(i) - key = value.GetName() - self.arrays.update({key: value}) - self.fieldtype = None - - - def hasArray(self, name): - """ - Check if the object contains an array with the requested name - - Parameters - ---------- - name : str - Name of the cell/point data array - - Returns - ------- - bool : True if the array exists. False otherwise - """ - if name in self.arrays.keys(): - return True - else: - return False - - - def getArray(self, name, sorted=False): - """ - Return the requested array from the cell data - - Parameters - ---------- - name : str - Name of the cell/point data array - sorted : bool - Sort with global ids \ - Require the presence of a global ids array - - Returns - ------- - numpy array : requested array - """ - array = None - - if self.hasArray(name): - array = VN.vtk_to_numpy(self.arrays[name]) - - if sorted: - ftype = self.fieldtype.split("Data")[0] - if self.hasArray(f"Global{ftype}Ids"): - gids = self.getCopyArray(f"Global{ftype}Ids") - array = array[np.argsort(gids)] - - return array - - - def getCopyArray(self, name, **kwargs): - """ - Return a copy of the requested array from the cell data - - Parameters - ---------- - name : str - Name of the cell/point data array - - Returns - ------- - numpy array : copy of the requested array - """ - - array = self.getArray(name, **kwargs) - - if array is not None: - array = copy.deepcopy(array) - - return array - - - def getVtkArray(self, name): - """ - Return the vtkDataArray requested - - Parameters - ---------- - name : str - Name of the cell/point data array - - Returns - ------- - numpy array : copy of the requested array - """ - if self.hasArray(name): - return self.arrays[name] - else: - return None - - - def setArray(self, name, value, overwrite=False): - """ - Return a copy of the requested array from the cell data - - Parameters - ---------- - name : str - Name of the cell data array - - Returns - ------- - numpy array : copy of the requested array - """ - if self.hasArray(name) and overwrite == False: - print(f"Warning! \n This dataset already contains a cell data array named {name}. Set the 'overwrite' parameter to True to bypass this warning") - else: - array = VN.vtk_to_numpy(self.arrays[name]) - array[:] = value[:] - - - -class VTKCellSpecifications(VTKFieldSpecifications): - """ - Contains the cell data information from a VTK Mesh - Inherits from VTKFieldSpecifications - """ - def __init__(self, celldata): - """ - Parameters - ---------- - celldata : vtk.vtkCellData - Cell data of the mesh - """ - assert(celldata.IsA("vtkCellData")) - super().__init__(fieldData=celldata) - self.fieldtype = "CellData" - - -class VTKPointSpecifications(VTKFieldSpecifications): - """ - Contains the point data information from a VTK Mesh - Inherits from VTKFieldSpecifications - - Parameters - ---------- - pointdata : vtk.vtkPointData - Point data of the mesh - - Attributes - --------- - arrays : dict - Dict containing the {name, vtkDataArray} of each point data array - """ - def __init__(self, pointdata): - """ - Parameters - ---------- - pointdata : vtk.vtkPointData - Point data of the mesh - """ - assert(pointdata.IsA("vtkPointData")) - super().__init__(fieldData=pointdata) - self.fieldtype = "PointData" diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py deleted file mode 100644 index f89a1e54..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py +++ /dev/null @@ -1,632 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import os -import vtk -from vtk.util import numpy_support as VN - -from .VtkFieldSpecifications import VTKCellSpecifications, VTKPointSpecifications -from ..model.utils.vtkUtils import cGlobalIds - - - -class VTKMesh: - """ - VTK format Mesh. Now handling .vtu .vts .pvts .pvtu - - Attributes - ---------- - meshfile : str - Mesh filename - vtktype : str - Format of the VTK mesh - bounds : tuple of int - Real bounds of the mesh (xmin, xmax, ymin, ymax, zmin, zmax) - numberOfPoints : int - Total number of points of the mesh - numberOfCells : int - Total number of cells of the mesh - isSet : bool - Whether or not the mesh properties have been set - hasLocator : bool - Whether or not the mesh cell locator has been initialized - """ - def __init__(self, meshfile): - """ - Parameters - ---------- - meshfile : str - Mesh filename - """ - self.meshfile = meshfile - self.vtktype = os.path.splitext(self.meshfile)[-1][1:] - - self.bounds = None - self.numberOfPoints = None - self.numberOfCells = None - - self.isSet = False - self.hasLocator = False - - - def getReader(self): - """Return the appropriate reader given the VTK format - - Returns - -------- - vtk.vtkXMLReader - Appropriate reader given the format of the VTK mesh file. - """ - - if self.vtktype == "vtu": - return vtk.vtkXMLUnstructuredGridReader() - elif self.vtktype == "vts": - return vtk.vtkXMLStructuredGridReader() - elif self.vtktype == "pvtu": - return vtk.vtkXMLPUnstructuredGridReader() - elif self.vtktype == "pvts": - return vtk.vtkXMLPStructuredGridReader() - else: - print("This VTK format is not handled.") - return None - - - def read(self): - """Read information from the VTK file - - Returns - -------- - vtk.vtkDataObject - General representation of VTK mesh data - """ - reader = self.getReader() - reader.SetFileName(self.meshfile) - reader.Update() - - return reader.GetOutput() - - - def setMeshProperties(self): - """Read and set as attributes the bounds, number of points and cells""" - data = self.read() - - self.bounds = data.GetBounds() - self.numberOfPoints = data.GetNumberOfPoints() - self.numberOfCells = data.GetNumberOfCells() - - self.isSet = True - - - def getBounds(self): - """ - Get the bounds of the mesh in the format: - (xmin, xmax, ymin, ymax, zmin, zmax) - - Returns - ------- - tuple or None - Bounds of the mesh - """ - return self.bounds - - - def getNumberOfPoints(self): - """ - Get the total number of points of the mesh - - Returns - ------- - int - Number of points - """ - return self.numberOfPoints - - - def getNumberOfCells(self): - """ - Get the total number of cells of the mesh - - Returns - ------- - int - Number of cells - """ - return self.numberOfCells - - - def getCellData(self): - """Read the cell data - - Returns - -------- - VTKCellSpecifications - Cell data information - """ - data = self.read() - - return VTKCellSpecifications(data.GetCellData()) - - - def getPointData(self): - """Read the point data - - Returns - -------- - VTKPointSpecifications - Point data information - """ - data = self.read() - - return VTKPointSpecifications(data.GetPointData()) - - - def getArray(self, name, dtype="cell", copy=False, sorted=False): - """ - Return a cell or point data array. If the file is a pvtu, the array is sorted with global ids - - Parameters - ----------- - name : str - Name of the vtk cell/point data array - dtype : str - Type of vtk data \ - `cell` or `point` - copy : bool - Return a copy of the requested array - Default is False - sorted : bool - Return the array sorted with respect to GlobalPointIds or GlobalCellIds - Default is False - - Returns - -------- - array : numpy array - Requested array - """ - assert dtype.lower() in ("cell", "point") - - if dtype.lower() == "cell": - fdata = self.getCellData() - else: - fdata = self.getPointData() - - if copy: - array = fdata.getCopyArray(name, sorted=sorted) - else: - array = fdata.getArray(name, sorted=sorted) - - return array - - - def extractMesh(self, center, srootname, dist=[None, None, None], comm=None, export=True): - """ - Extract a rectangular submesh such that for each axis we have the subax: [center-dist, center+dist] - - Parameters - --------- - center : 3d float - Requested center of the subbox - srootname : str - Submesh root filename - dist : 3d float - Distance to the center in each direction - comm : MPI.COMM_WORLD - MPI communicator - - Returns - ------- - VTKMesh - Submesh extracted - """ - assert self.vtktype in ("vtu", "pvtu", "vts", "pvts") - vtktype = self.vtktype[-3:] - sfilename = ".".join((srootname, vtktype)) - - if comm is None or comm.Get_rank() == 0: - if not self.isSet: - self.setMeshProperties() - minpos = [] - maxpos = [] - - for i in range(3): - xmin, d = self.getSubAx(center[i], dist[i], ax=i+1) - minpos.append(xmin) - maxpos.append(xmin+d) - - data = self.read() - submesh = VTKSubMesh(sfilename, data, minpos, maxpos, create=export) - - else: - submesh = VTKMesh(sfilename) - - # if file creation, wait for rank 0 to finish - if export: - info = "Done" - comm.bcast(info, root=0) - - return submesh - - - def getSubAx(self, center, dist, ax): - """ - Return the min and max positions in the mesh given the center, distance and ax considered. If the 2*distance if greater than the bounds, the min/max is the corresponding mesh bound. - - Parameters - ---------- - center : float - Central position considered - dist : float - Max distance requested - ax : int - Ax to consider (1, 2, 3) - - Returns - ------- - min, max : float - Min and Max positions - """ - assert(type(ax) == int) - - bounds = [self.bounds[(ax-1)*2], self.bounds[ax*2-1]] - - if dist is not None: - dist = abs(dist) - ox = max(bounds[0], center-dist) - x = min(bounds[1]-ox, 2*dist) - else: - ox = bounds[0] - x = bounds[1] - - return ox, x - - - def getNumberOfBlocks(self): - """Return the number of blocks of a mesh.""" - if self.vtktype in ["pvtu", "pvts"]: - with open(self.meshfile) as ff: - nb = 0 - for line in ff: - m = line.split() - if m[0] == ' high are set to high - comm : MPI_COMM - MPI communicators - """ - root = 0 - rank = comm.Get_rank() - size = comm.Get_size() - - x = None - if rank == root: - with h5py.File( filename, 'r' ) as f: - x = f["velocity"][:] - - imin = np.where( x < low )[0] - imax = np.where( x > high )[0] - x[imin] = low - x[imax] = high - - if self.model == "1/c2": - x = np.sqrt(1/x) - elif self.model == "1/c": - x = 1/x - elif self.model == "c": - pass - else: - raise ValueError("Not implemented") - - startModel = self.bcastField( x, comm ) - - self.updateVelocityModel( startModel ) - - - def updateVelocityModel(self, vel, **kwargs): - """ - Update velocity value in GEOS - - Parameters - ---------- - vel : array - Values for velocity field - """ - super().updateVelocityModel(vel, velocityName="acousticVelocity", **kwargs) - - - def setConstantVelocityModel(self, vel, velFieldName="acousticVelocity", **kwargs): - """ - Update velocity value in GEOS, using a constant value. - - Parameters - ---------- - vel : float - Value for velocity field - """ - prefix = self._getPrefixPath(**kwargs) - - velocity = self.solver.get_wrapper(prefix+velFieldName).value() - velocity.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - velocity.to_numpy()[:] = vel - - - def getVelocityModel(self, filterGhost=False, **kwargs): - """ - Get the velocity values - - Parameters - ----------- - filterGhost : bool - Filter the ghost ranks - - Returns - ------- - velocity : numpy array - Velocity values - """ - velocity = super().getVelocityModel(velocityName="acousticVelocity", filterGhost=filterGhost, **kwargs) - - return velocity - - - def updateDensityModel(self, density, **kwargs): - """ - Update density values in GEOS - - Parameters - ----------- - density : array - New values for the density - """ - super().updateDensityModel( density=density, densityName="acousticDensity", **kwargs ) - - - def getGradient(self, filterGhost=False, **kwargs): - """Get the local rank gradient value - - Returns - -------- - partialGradient : Numpy Array - Array containing the element id list for the local rank - """ - partialGradient = self.getField("partialGradient", **kwargs) - - if filterGhost: - partialGradient = self.filterGhostRank(partialGradient, **kwargs) - - return partialGradient - - - def computePartialGradient( self, shotId, minDepth, comm, gradDirectory="partialGradient", **kwargs ): - """Compute the partial Gradient - - Parameters - ----------- - shotId : string - Number of the shot as string - minDepth : float - Depth at which gradient values are kept, otherwise it is set to 0. - NOTE : this is done in this routine to avoid storage \ - of elementCenter coordinates in the .hdf5 \ - but might be problem for WolfeConditions later on \ - if minDepth is too large - comm : MPI_COMM - MPI communicators - gradDirectory : str, optional - Partial gradient directory \ - Default is `partialGradient` - """ - rank = comm.Get_rank() - size = comm.Get_size() - root = 0 - - # Get local gradient - grad = self.getGradient( **kwargs ) - if self.model == "1/c2": - x = self.getVelocityModel( **kwargs ) - grad = -( x * x * x / 2 ) * grad - elif self.model == "1/c": - x = self.getVelocityModel( filterGhost=True, **kwargs ) - grad = - x * x * grad - elif self.model == "c": - pass - else: - raise ValueError("Not implemented") - - grad = grad.astype( np.float64 ) - - zind = np.where(self.getElementCenter(**kwargs)[:,2] < minDepth)[0] - grad[zind] = 0.0 - - # Gather global gradient - gradFull, ntot = self.gatherField( field=grad, comm=comm, root=root ) - - if rank == root: - os.makedirs(gradDirectory, exist_ok=True) - - with h5py.File(f"{gradDirectory}/partialGradient_"+shotId+".hdf5", 'w') as h5p: - - h5p.create_dataset("partialGradient", data=np.zeros(ntot), chunks=True, maxshape=(ntot,)) - h5p["partialGradient"][:] = self.dtWaveField * gradFull - - shutil.move(f"{gradDirectory}/partialGradient_"+shotId+".hdf5", f"{gradDirectory}/partialGradient_ready_"+shotId+".hdf5") - - comm.Barrier() - - - def getPressureAtReceivers(self): - """ - Get the pressure values at receivers coordinates - - Returns - ------ - numpy Array : Array containing the pressure values at all time step at all receivers coordinates - """ - pressureAtReceivers = self.solver.get_wrapper("pressureNp1AtReceivers").value() - - return pressureAtReceivers.to_numpy() - - - def getFullPressureAtReceivers(self, comm): - """Return all pressures at receivers values on all ranks - Note that for a too large 2d array this may not work. - - Parameters: - ----------- - comm : MPI_COMM - MPI communicators - """ - rank = comm.Get_rank() - - allPressure = comm.gather(self.getPressureAtReceivers(), root=0) - pressure = np.zeros(self.getPressureAtReceivers().shape) - - if rank == 0: - for p in allPressure: - for i in range(p.shape[1]): - if any(p[1:, i])==True: - pressure[:, i] = p[:, i] - - pressure = comm.bcast(pressure, root=0) - - return pressure - - - def resetWaveField(self, **kwargs): - """Reinitialize all pressure values on the Wavefield to zero in GEOSX""" - - self.geosx.get_wrapper("Solvers/"+self.name+"/indexSeismoTrace").value()[0] = 0 - meshName = self._getMeshName() - discretization = self._getDiscretization() - - nodeManagerPath = f"domain/MeshBodies/{meshName}/meshLevels/{discretization}/nodeManager/" - - - if self.type == "AcousticSEM": - for ts in ("nm1", "n", "np1"): - pressure = self.geosx.get_wrapper(nodeManagerPath + f"pressure_{ts}").value() - pressure.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - pressure.to_numpy()[:] = 0.0 - - elif self.type == "AcousticFirstOrderSEM": - pressure_np1 = self.geosx.get_wrapper(nodeManagerPath + "pressure_np1").value() - pressure_np1.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - pressure_np1.to_numpy()[:] = 0.0 - - prefix = self._getPrefixPath(**kwargs) - for component in ("x", "y", "z"): - velocity = self.geosx.get_wrapper(prefix + f"velocity_{component}").value() - velocity.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - velocity.to_numpy()[:] = 0.0 - - - def resetPressureAtReceivers(self): - """Reinitialize pressure values at receivers to 0 - """ - pressure = self.solver.get_wrapper("pressureNp1AtReceivers").value() - pressure.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - pressure.to_numpy()[:] = 0.0 - - - def getWaveField(self): - return self.getPressureAtReceivers()[:,:-1] - - - def getFullWaveFieldAtReceivers(self, comm): - return self.getFullPressureAtReceivers(comm)[:,:-1] diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py deleted file mode 100644 index 34acbc77..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py +++ /dev/null @@ -1,319 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import numpy as np -import ctypes as ct - -import pygeosx - -from .WaveSolver import WaveSolver - - -class ElasticSolver(WaveSolver): - """ - ElasticSolver Object containing all methods to run ElasticSEM simulation with GEOSX - - Attributes - ----------- - dt : float - Time step for simulation - minTime : float - Min time to consider - maxTime : float - End time to consider - dtSeismo : float - Time step to save pressure for seismic trace - minTimeSim : float - Starting time of simulation - maxTimeSim : float - End Time of simulation - dtWaveField : float - Time step to save fields - sourceType : str - Type of source - sourceFreq : float - Frequency of the source - name : str - Solver name - type : str - Type of solver - Default is None - geosxArgs : GeosxArgs - Object containing GEOSX launching options - geosx : pygeosx instance - - alreadyInitialized : bool - Flag indicating if the initial conditions have been applied yet - firstGeosxInit : bool - Flag for initialize or reinitialize - collectionTargets : list - Output targets for geosx - hdf5Targets : list - HDF5 output targets for geosx - vtkTargets : list - VTK output targets for geosx - """ - def __init__(self, - dt=None, - minTime=0, - maxTime=None, - dtSeismo=None, - dtWaveField=None, - sourceType=None, - sourceFreq=None, - **kwargs): - """ - Parameters - ---------- - dt : float - Time step for simulation - minTime : float - Starting time of simulation - Default is 0 - maxTime : float - End Time of simulation - dtSeismo : float - Time step to save pressure for seismic trace - dtWaveField : float - Time step to save fields - sourceType : str - Type of source - Default is None - sourceFreq : float - Frequency of the source - Default is None - kwargs : keyword args - geosx_argv : list - GEOSX arguments or command line as a splitted line - """ - super().__init__(dt=dt, - minTime=minTime, - maxTime=maxTime, - dtSeismo=dtSeismo, - dtWaveField=dtWaveField, - sourceType=sourceType, - sourceFreq=sourceFreq, - **kwargs) - - def initialize(self, rank=0, xml=None): - super().initialize(rank, xml) - try: - useDAS = self.xml.getAttribute(parentElement=self._getType(), attributeTag="useDAS") - - except AttributeError: - useDAS = None - - if useDAS == "none": - try: - linearGEO = bool(self.xml.getAttribute(self._getType(), "linearDASGeometry")) - except AttributeError: - linearGEO = False - - if linearGEO is True: - self.useDAS = True - - - def __repr__(self): - string_list = [] - string_list.append("Solver type : " + type(self).__name__ + "\n") - string_list.append("dt : " + str(self.dt) + "\n") - string_list.append("maxTime : " + str(self.maxTime) + "\n") - string_list.append("dtSeismo : " + str(self.dtSeismo) + "\n") - string_list.append("Outputs : " + str(self.hdf5Targets) + "\n" + str(self.vtkTargets) + "\n") - rep = "" - for string in string_list: - rep += string - - return rep - - - def updateVelocityModel(self, vel, component, **kwargs): - """ - Update velocity value in GEOS - - Parameters - ---------- - vel : float/array - Value(s) for velocity field - component : str - Vs or Vp - """ - assert component.lower() in ("vs", "vp"), "Only Vs or Vp component accepted" - super().updateVelocityModel(vel, velocityName="elasticVelocity"+component.title(), **kwargs) - - - def getVelocityModel(self, component, filterGhost=False, **kwargs): - """ - Get the velocity values - - Parameters - ----------- - component : str - Vs or Vp - filterGhost : bool - Filter the ghost ranks - - Returns - ------- - velocity : numpy array - Array containing the velocity values - """ - assert component.lower() in ("vs", "vp"), "Only Vs or Vp component accepted" - - velocity = super().getVelocityModel(velocityName="elasticVelocity"+component.title(), filterGhost=filterGhost, **kwargs) - - return velocity - - - def getDensityModel(self, filterGhost=False, **kwargs): - """ - Get the density values - - Parameters - ----------- - filterGhost : bool - Filter the ghost ranks - - Returns - -------- - density : numpy array - Array containing the density values - """ - density = self.getField("elasticDensity", filterGhost=filterGhost, **kwargs) - - return density - - - def updateDensityModel(self, density, **kwargs): - """ - Update density values in GEOS - - Parameters - ----------- - density : array - New values for the density - """ - super().updateDensityModel( density=density, densityName="elasticDensity", **kwargs ) - - - def getDASSignalAtReceivers(self): - """ - Get the DAS signal values at receivers coordinates - - Returns - -------- - dassignal : numpy array - Array containing the DAS signal values at all time step at all receivers coordinates - """ - if self.type != "ElasticSEM": - raise TypeError(f"DAS signal not implemented for solver of type {self.type}.") - else: - dassignal = self.solver.get_wrapper(f"dasSignalNp1AtReceivers").value().to_numpy() - - return dassignal - - - def getDisplacementAtReceivers(self, component="X"): - """ - Get the displacement values at receivers coordinates for a given direction - - Returns - -------- - displacement : numpy array - Array containing the displacements values at all time step at all receivers coordinates - """ - assert component.upper() in ("X", "Y", "Z") - if self.type == "ElasticFirstOrderSEM": - displacement = self.solver.get_wrapper(f"displacement{component.lower()}Np1AtReceivers").value().to_numpy() - elif self.type == "ElasticSEM": - displacement = self.solver.get_wrapper(f"displacement{component.upper()}Np1AtReceivers").value().to_numpy() - - return displacement - - - def getAllDisplacementAtReceivers(self): - """ - Get the displacement for the x, y and z directions at all time step and all receivers coordinates - - Returns - -------- - displacementX : numpy array - Component X of the displacement - displacementY : numpy array - Component Y of the displacement - displacementZ : numpy array - Component Z of the displacement - """ - displacementX = self.getDisplacementAtReceivers("X") - displacementY = self.getDisplacementAtReceivers("Y") - displacementZ = self.getDisplacementAtReceivers("Z") - - return displacementX, displacementY, displacementZ - - - def resetWaveField(self, **kwargs): - """Reinitialize all displacement values on the Wavefield to zero in GEOSX""" - - self.geosx.get_wrapper("Solvers/"+self.name+"/indexSeismoTrace").value()[0] = 0 - meshName = self._getMeshName() - discretization = self._getDiscretization() - - nodeManagerPath = f"domain/MeshBodies/{meshName}/meshLevels/{discretization}/nodeManager/" - - if self.type == "ElasticSEM": - for component in ("x", "y", "z"): - for ts in ("nm1", "n", "np1"): - displacement = self.geosx.get_wrapper(nodeManagerPath+f"displacement{component}_{ts}").value() - displacement.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - displacement.to_numpy()[:] = 0.0 - - elif self.type == "ElasticFirstOrderSEM": - component = ("x", "y", "z") - for c in component: - displacement_np1 = self.geosx.get_wrapper(nodeManagerPath+f"displacement{c}_np1").value() - displacement_np1.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - displacement_np1.to_numpy()[:] = 0.0 - - prefix = self._getPrefixPath(**kwargs) - for i, c in enumerate(component): - for j in range(i, len(component)): - cc = c + component[j] - - sigma = self.solver.get_wrapper(prefix+f"stresstensor{cc}").value() - sigma.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - sigma.to_numpy()[:] = 0.0 - - - def resetDisplacementAtReceivers(self): - """Reinitialize displacement values at receivers to 0 - """ - for component in ("X", "Y", "Z"): - displacement = self.solver.get_wrapper(f"displacement{component}Np1AtReceivers").value() - displacement.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - displacement.to_numpy()[:] = 0.0 - - - def getWaveField( self ): - if self.useDAS: - return self.getDASSignalAtReceivers() - else: - return self.getAllDisplacementAtReceivers() - - - def getFullWaveFieldAtReceivers(self, comm): - print( "This method is not implemented yet" ) diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py deleted file mode 100644 index ebf16fea..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py +++ /dev/null @@ -1,104 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -from .Solver import Solver - - -class GeomechanicsSolver(Solver): - """ - Geomechanics solver object containing methods to run geomechanics simulations in GEOS - - Attributes - ----------- - xml : XML - XML object containing parameters for GEOSX initialization - initialDt : float - Time step for the simulation - maxTime : float - End time of simulation - name : str - Solver name - type : str - Type of solver - Default is None - geosxArgs : GeosxArgs - Object containing GEOSX launching options - geosx : pygeosx.Group - pygeosx initialize instance - alreadyInitialized : bool - Flag indicating if the initial conditions have been applied yet - firstGeosxInit : bool - Flag for initialize or reinitialize - collectionTargets : list - Output targets for geosx - hdf5Targets : list - HDF5 output targets for geosx - vtkTargets : list - VTK output targets for geosx - """ - def __init__(self, - dt=None, - maxTime=None, - **kwargs): - """ - Parameters - ----------- - initialDt : float - Initial time step \ - Default is None - maxTime : float - End time of simulation \ - Default is None - """ - super().__init__(**kwargs) - - self.dt = dt - self.maxTime = maxTime - - - def _setTimeVariables(self): - """Set the time variables attributes""" - # Set up variables from xml - events = self.xml.events - - if self.maxTime is None: - self.maxTime = float(events["maxTime"]) - else: - self.updateTimeVariable("maxTime") - - if self.dt is None: - for event in events["PeriodicEvent"]: - if isinstance(event, dict): - poroName = "poromechanicsSolver" - if event["target"] == "/Solvers/" + poroName: - self.dt = float(event['forceDt']) - else: - if event == "target" and events["PeriodicEvent"]["target"] == "/Solvers/" + self.name: - self.dt = float(events["PeriodicEvent"]["forceDt"]) - else: - self.updateTimeVariable("dt") - - - def updateTimeVariable(self, variable): - """Overwrite a time variable in GEOS""" - if variable == "maxTime": - assert hasattr(self, "maxTime") - self.geosx.get_wrapper("Events/maxTime").value()[0] = self.maxTime - elif variable == "dt": - assert hasattr(self, "dt") - self.geosx.get_wrapper("Events/solverApplications/forceDt").value()[0] - - - def initialize( self, rank=0, xml=None ): - super().initialize( rank, xml, stype="SinglePhasePoromechanics" ) diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py deleted file mode 100644 index 48f7a438..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py +++ /dev/null @@ -1,216 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import numpy as np - -from .Solver import Solver - - -class ReservoirSolver(Solver): - """ - Reservoir solver object containing methods to run reservoir simulations in GEOS - - Attributes - ----------- - xml : XML - XML object containing parameters for GEOSX initialization - initialDt : float - Time step for the simulation - maxTime : float - End time of simulation - name : str - Solver name - type : str - Type of solver - Default is None - geosxArgs : GeosxArgs - Object containing GEOSX launching options - geosx : pygeosx.Group - pygeosx initialize instance - alreadyInitialized : bool - Flag indicating if the initial conditions have been applied yet - firstGeosxInit : bool - Flag for initialize or reinitialize - collectionTargets : list - Output targets for geosx - hdf5Targets : list - HDF5 output targets for geosx - vtkTargets : list - VTK output targets for geosx - """ - def __init__(self, - initialDt=None, - maxTime=None, - **kwargs): - """ - Parameters - ----------- - initialDt : float - Initial time step \ - Default is None - maxTime : float - End time of simulation \ - Default is None - """ - super().__init__(**kwargs) - - self.initialDt = initialDt - self.maxTime = maxTime - - self.isCoupled = kwargs.get("coupled", False) - self.dtSeismic4D = None - - - def _setTimeVariables(self): - """Set the time variables attributes""" - # Set up variables from xml - events = self.xml.events - outputs = self.xml.outputs - - if self.isCoupled: - if self.dtSeismic4D is None: - for event in events["PeriodicEvent"]: - if event["name"] == "seismic4D": - self.dtSeismic4D = float(event["timeFrequency"]) - else: - self.updateTimeVariable("dtSeismic4D") - - if self.maxTime is None: - self.maxTime = float(events["maxTime"]) - else: - self.updateTimeVariable("maxTime") - - fvm = self.xml.solvers[self.type] - if self.initialDt is None: - for k, v in fvm.items(): - if k == "initialDt": - self.initialDt = np.array(v, dtype=float) - else: - self.updateTimeVariable("initialDt") - - - def updateTimeVariable(self, variable): - """Overwrite a time variable in GEOS""" - if variable == "maxTime": - assert hasattr(self, "maxTime") - self.geosx.get_wrapper("Events/maxTime").value()[0] = self.maxTime - elif variable == "initialDt": - assert hasattr(self, "initialDt") - self.geosx.get_wrapper(f"/Solvers/{self.name}/initialDt").value()[0] = self.initialDt - elif variable == "dtSeismic4D": - assert hasattr(self, "dtSeismic4D") - self.geosx.get_wrapper("Events/seismic4D/timeFrequency").value()[0] = self.dtSeismic4D - - - def execute(self, time): - """ - Do one solver iteration - - Parameters - ---------- - time : float - Current time of simulation - dt : float - Timestep - """ - self.solver.execute(time, self.initialDt) - - - def getTimeStep(self): - """ - Get the value of `initialDt` variable from GEOS - - Returns - -------- - float - initialDt value - """ - return self.solver.get_wrapper("initialDt").value()[0] - - - def updateTimeStep(self): - """Update the attribute value with the one from GEOS""" - self.initialDt = self.getTimeStep() - - - def getPressure(self, **kwargs): - """ - Get the local pressure - - Returns - -------- - pressure : numpy array - Local pressure - """ - pressure = self.getField("pressure", **kwargs) - - return pressure - - - def getDeltaPressure(self, **kwargs): - """ - Get the local delta pressure - - Returns - -------- - deltaPressure : numpy array - Local delta pressure - """ - deltaPressure = self.getField("deltaPressure", **kwargs) - - return deltaPressure - - - def getPhaseVolumeFractionGas(self, **kwargs): - """ - Get the local gas phase volume fraction - - Returns - -------- - phaseVolumeFractionGas : numpy array - Local gas phase volume fraction - """ - phaseVolumeFraction = self.getField("phaseVolumeFraction", **kwargs) - phaseVolumeFractionGas = np.ascontiguousarray(phaseVolumeFraction[:,0]) - - return phaseVolumeFractionGas - - - def getPhaseVolumeFractionWater(self, **kwargs): - """ - Get the local water phase volume fraction - - Returns - -------- - phaseVolumeFractionWater : numpy array - Local water phase volume fraction - """ - phaseVolumeFraction = self.getField("phaseVolumeFraction", **kwargs) - phaseVolumeFractionWater = np.ascontiguousarray(phaseVolumeFraction[:,1]) - - return phaseVolumeFractionWater - - - def getRockPorosity(self, **kwargs): - """ - Get the local rock porosity - - Returns - -------- - rockPorosity : numpy array - Local rock porosity - """ - rockPorosity = self.getField("rockPorosity_referencePorosity", **kwargs) - - return rockPorosity diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py deleted file mode 100644 index a2afd4d9..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py +++ /dev/null @@ -1,697 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -import os -import sys -import numpy as np - -import pygeosx -from ..input.Xml import XML -from ..input.GeosxArgs import GeosxArgs - -from mpi4py import MPI - -class Solver: - """ - Solver class containing the main methods of a GEOS solver - """ - - def __init__(self, **kwargs): - self.alreadyInitialized = False - self.type = None - - argv = kwargs.get("geosx_argv", sys.argv) - self.geosxArgs = GeosxArgs(argv) - self.xml = None - - - def initialize(self, rank=0, xml=None, **kwargs): - """Initialization or reinitialization of GEOSX - - Parameters - ---------- - rank : int - Process rank - xml : XML - XML object containing parameters for GEOSX initialization. - Only required if not set in the __init__ OR if different from it - """ - if xml: - self.updateXml(xml) - else: - if self.xml is None: - try: - self.xml = XML(self.geosxArgs.options["xml"]) - except: - raise ValueError("You need to provide a xml input file") - - - if not self.alreadyInitialized: - # if GEOS is UNINITIALIZED - if self.getGEOSState() == 0: - self.geosx = pygeosx.initialize(rank, self.geosxArgs.getCommandLine()) - self.alreadyInitialized = True - - # elif GEOS is INITIALIZED OR READY TO RUN - elif self.getGEOSState() in (1, 2): - self.geosx = pygeosx.reinit(self.geosxArgs.getCommandLine()) - self.alreadyInitialized = True - - # else if COMPLETED state - else: - print(f"The current state of GEOS does not allow for initialization\nCurrent state: {self.getGEOSState()}") - - self.name = self._getName() - stype = kwargs.get( "stype", None ) - self._setSolverGroup( stype ) - - # Check all time variables are set/ set them from xml - self._setTimeVariables() - - # Set the outputs collections/targets - self.__setOutputs() - - - def _setSolverGroup( self, stype=None ): - if stype is None: - if self.name is not None: - if isinstance( self.name, str ): - self.solver = self.geosx.get_group( "/Solvers/" + self.name ) - else: - name = self._getName( stype=stype ) - self.solver = self.geosx.get_group( "/Solvers/" + name ) - - - def getGEOSState(self): - """ - Return the current GEOS state - - Returns - --------- - int - GEOS state - 0 : UNINITIALIZED - 1 : INITIALIZED - 2 : READY TO RUN - 3 : COMPLETED - """ - return pygeosx.getState() - - - def _setTimeVariables(self): - """Initialize the time variables. Specific to each solver""" - pass - - - def __setOutputs(self): - if hasattr(self.xml, "outputs"): - outputs = self.xml.outputs - self.__setOutputTargets(outputs) - - self.collections = [] - self.hdf5Outputs = [] - self.vtkOutputs = [] - - if not outputs == None: - # Set the collections - for target in self.collectionTargets: - self.collections.append(self.geosx.get_group(target)) - - for target in self.hdf5Targets: - self.hdf5Outputs.append(self.geosx.get_group(target)) - - for target in self.vtkTargets: - self.vtkOutputs.append(self.geosx.get_group(target)) - - - def __setOutputTargets(self, outputs): - """ - Set the list of outputs targets - - Parameters - ----------- - outputs : dict - Dictionnary extracted from the xml input file - """ - self.collectionTargets = [] - self.hdf5Targets = [] - self.vtkTargets = [] - - if outputs: - if isinstance(list(outputs.values())[0], list): - if "TimeHistory" in outputs.keys(): - for hdf5 in outputs["TimeHistory"]: - self.collectionTargets.append(hdf5['sources'].strip("{} ")) - self.hdf5Targets.append("Outputs/"+hdf5['name']) - - if "VTK" in outputs.keys(): - for vtk in outputs["VTK"]: - self.vtkTargets.append("Outputs/"+vtk['name']) - - else: - if "TimeHistory" in list(outputs.keys()): - hdf5 = outputs["TimeHistory"] - self.collectionTargets.append(hdf5['sources'].strip("{} ")) - self.hdf5Targets.append("Outputs/"+hdf5['name']) - - if "VTK" in list(outputs.keys()): - vtk = outputs["VTK"] - self.vtkTargets.append("Outputs/"+vtk['name']) - - - def _getType(self): - """ - Set the type of solver given in the xml 'Solvers' block - - Raises - ------- - ValueError : if no solver is provided in the XML - """ - if self.xml is not None: - typesOfSolvers = self.xml.getSolverType() - - if len( typesOfSolvers ) == 1: - self.type = typesOfSolvers[0] - elif len( typesOfSolvers ) > 1: - self.type = typesOfSolvers - else: - raise ValueError("You must provide a Solver in the XML input file") - - - def _getName(self, stype=None): - """ - Get the solver 'name' attribute from the xml - - Returns - ------- - str or list of str - Solver name from the xml \ - List of solvers name if several solvers in xml - """ - if self.xml is not None: - if stype is None: - if self.type is None: - self._getType() - stype = self.type - - # Check only one type of solver - if isinstance(stype, str): - return self.xml.solvers[stype]["name"] - - elif isinstance( stype, list ): - return [ self.xml.solvers[solvertype]["name"] for solvertype in stype ] - - - def _getMeshName(self): - """ - Get the mesh 'name' attribute from the xml - - Returns - ------- - str - Mesh name from the xml - """ - if self.xml is not None: - meshes = [m for m in self.xml.mesh] - if len(meshes) <= 1: - return self.xml.mesh[meshes[0]]["name"] - - - def _getDiscretization(self): - """Get discretization from the XML - - Returns - -------- - discretization : str - Name of the discretization method - """ - if self.xml is not None: - if self.type is None: - self._getType() - - if isinstance(self.type, str): - return self.xml.solvers[self.type]["discretization"] - elif isinstance(self.type, list): - return [self.xml.solvers[solverType]["discretization"] for solverType in self.type] - - - def _getTargetRegion(self): - """ - Get the target region name from the xml - - Returns - ------- - str - target region from the xml - """ - if self.xml is not None: - if self.type is None: - self._getType() - - if isinstance(self.type, str): - targetRegionRaw = self.xml.solvers[self.type]["targetRegions"] - targetRegion = targetRegionRaw.strip("{ }").split(",") - - if len(targetRegion) <= 1: - return targetRegion[0] - - - def _getCellBlock(self): - """ - Get the cell blocks names from the xml - - Returns - ------- - str - cell blocks from the xml - """ - if self.xml is not None: - cellElementRegion = self.xml.elementRegions["CellElementRegion"][0] - cellBlocks = cellElementRegion["cellBlocks"].strip("{ }").split(",") - - if len(cellBlocks) <= 1: - return cellBlocks[0] - - - def reinitSolver(self): - """Reinitialize Solver""" - self.solver.reinit() - - - def applyInitialConditions(self): - """Apply the initial conditions after GEOS (re)initialization""" - if self.getGEOSState() == 1: - pygeosx.apply_initial_conditions() - - - - def finalize(self): - """Terminate GEOSX""" - pygeosx._finalize() - - - def updateXml(self, xml): - """ - Update XML - - Parameters - ----------- - xml : XML - XML object corresponding to GEOSX input - """ - self.xml = xml - - if self.geosxArgs.updateArg("xml", xml.filename): - self.alreadyInitialized = False - - - def updateHdf5OutputsName(self, directory, filenames, reinit=False): - """ - Overwrite GEOSX hdf5 Outputs paths that have been read in the XML. - - Parameters - ---------- - list_of_output : list of str - List of requested output paths - reinit : bool - Perform reinitialization or not. Must be set to True if called after applyInitialConditions() - """ - - if not len(self.hdf5Outputs): - raise ValueError("No HDF5 Outputs specified in XML.") - else: - for i in range(len(filenames)): - os.makedirs(directory, exist_ok=True) - - self.hdf5Outputs[i].setOutputName(os.path.join(directory, filenames[i])) - if reinit: - self.hdf5Outputs[i].reinit() - - - def updateVtkOutputsName(self, directory): - """ - Overwrite GEOSX vtk Outputs paths that have been read in the XML. - - Parameters - ---------- - list_of_output : list of str - List of vtk output paths - reinit : bool - Perform reinitialization or not. Must be set to True if called after applyInitialConditions() - """ - if not len(self.vtkOutputs): - pass - else: - self.vtkOutputs[0].setOutputDir(directory) - - - def execute(self, time): - """ - Do one solver iteration - - Parameters - ---------- - time : float - Current time of simulation - """ - - self.solver.execute(time, self.dt) - - - def cleanup(self, time): - """ - Finalize simulation. Also triggers write of leftover seismogram data - - Parameters - ---------- - time : float - Current time of simulation - """ - self.solver.cleanup(time) - - - def outputVtk(self, time): - """ - Trigger the VTK output - - Parameters - ---------- - time : float - Current time of simulation - """ - for vtkOutput in self.vtkOutputs: - vtkOutput.output(time, self.dt) - - - def _getPrefixPath(self, targetRegion=None, meshName=None, cellBlock=None): - """ - Return the prefix path to get wrappers or fields in GEOS - - Parameters - ----------- - targetRegion : str, optional - Name of the target Region \ - Default value is taken from the xml - meshName : str, optional - Name of the mesh \ - Default value is taken from the xml - cellBlock : str, optional - Name of the cell blocks \ - Default value is taken from the xml - - Returns - ------- - prefix : str - Prefix path - - Raises - ------- - AssertionError : if the variables 'targetRegion', 'meshName' \ - or `cellBlock` have multiple or no values - """ - if targetRegion is None: - targetRegion = self._getTargetRegion() - if meshName is None: - meshName = self._getMeshName() - if cellBlock is None: - cellBlock = self._getCellBlock() - - discretization=self._getDiscretization() - if discretization is None: - discretization="Level0" - assert None not in (targetRegion, meshName, cellBlock, discretization), "No values or multiple values found for `targetRegion`, `meshName` and `cellBlock` arguments" - - prefix = os.path.join("/domain/MeshBodies", meshName, "meshLevels", discretization, "ElementRegions/elementRegionsGroup", targetRegion, "elementSubRegions", cellBlock, "") - return prefix - - - def getField(self, fieldName, **kwargs): - """ - Get the requested field as numpy array - - Parameters - ----------- - fieldName : str - Name of the field in GEOSX - - Returns - ------- - field : np.array - Field requested - """ - prefix = self._getPrefixPath(**kwargs) - field = self.solver.get_wrapper(prefix+fieldName).value() - - return field.to_numpy() - - - def getElementCenter(self, filterGhost=False, **kwargs): - """ - Get element center position as numpy array - - Returns - ------- - elementCenter : array-like - Element center coordinates - """ - elementCenter = self.getField("elementCenter", **kwargs) - - if filterGhost: - elementCenter = self.filterGhostRank(elementCenter, **kwargs) - - return elementCenter - - - def getElementCenterZ(self, **kwargs): - """ - Get the z coordinate of the element center - - Returns - ------- - elementCenterZ : array-like - Element center z coordinates - """ - elementCenter = self.getField("elementCenter", **kwargs) - elementCenterZ = np.ascontiguousarray(elementCenter[:,2]) - - return elementCenterZ - - - def getGhostRank(self, **kwargs): - """ - Get the local ghost ranks - - Returns - ------- - ghostRank : array-like - Local ghost ranks - """ - ghostRank = self.getField("ghostRank", **kwargs) - - return ghostRank - - - def getLocalToGlobalMap(self, filterGhost=False, **kwargs): - """ - Get the local rank element id list - - Returns - ------- - Numpy Array : Array containing the element id list for the local rank - """ - localToGlobalMap = self.getField("localToGlobalMap", **kwargs) - - if filterGhost: - localToGlobalMap = self.filterGhostRank(localToGlobalMap, **kwargs) - - return localToGlobalMap - - - def gatherField(self, field, comm, root=0, **kwargs): - """ - Gather a full GEOS field from all local ranks - - Parameters - ----------- - field : numpy array - Local field - comm : MPI.COMM_WORLD - MPI communicator - root : int - MPI rank used for the gather \ - Default is rank 0 - """ - assert isinstance(root, int) - assert root < comm.Get_size() - - rank = comm.Get_rank() - - ghostRank = self.getGhostRank(**kwargs) - localToGlobalMap = self.getLocalToGlobalMap(**kwargs) - - # Prepare buffer - nlocalElements = ghostRank.shape[0] - nmax = np.zeros( 1 ) - nmax[0] = np.max( localToGlobalMap ) # max global number of elements - - comm.Barrier() - comm.Allreduce( MPI.IN_PLACE, nmax, op=MPI.MAX ) - ntot = round( nmax[0] + 1 ) - - if rank != root: - fullField = None - nrcv = nlocalElements - comm.send(nrcv, dest=root, tag=1) - comm.Send(field, dest=root, tag=2) - comm.Send(ghostRank, dest=root, tag=3) - comm.Send(localToGlobalMap, dest=root, tag=4) - - else: - fullField = np.full( (ntot), fill_value=np.nan ) - jj = np.where( ghostRank < 0 )[0] - fullField[localToGlobalMap[jj]] = field[jj] - - for r in range( comm.Get_size() ): - if r != root: - nrcv = comm.recv(source=r, tag=1) - - fieldRcv = np.zeros(nrcv, dtype=np.float64) - ghostRankRcv = np.zeros(nrcv, dtype=np.int32) - localToGlobalMapRcv = np.zeros(nrcv, dtype=np.int64) - - comm.Recv(fieldRcv, source=r, tag=2) - comm.Recv(ghostRankRcv, source=r, tag=3) - comm.Recv(localToGlobalMapRcv, source=r, tag=4) - - jj = np.where ( ghostRankRcv < 0 )[0] - - fullField[ localToGlobalMapRcv[jj]] = fieldRcv[jj] - comm.Barrier() - return fullField, ntot - - - def bcastField( self, fullField, comm, root=0, **kwargs ): - """ - Broadcast a field to local ranks with GEOS local to global map - - Parameters - ----------- - fullField : numpy array - Full field - comm : MPI.COMM_WORLD - MPI communicator - root : int - MPI rank used for the gather \ - Default is rank 0 - - Returns - -------- - field : numpy array - Local field - """ - rank = comm.Get_rank() - size = comm.Get_size() - - ghostRank = self.getGhostRank( **kwargs ) - localToGlobalMap = self.getLocalToGlobalMap( **kwargs ) - nlocalElements = ghostRank.shape[0] - - field = np.zeros( nlocalElements ) - - if rank == root: - jj = np.where(ghostRank < 0)[0] - field[jj] = fullField[localToGlobalMap[jj]] - - for r in range( size ): - if r != root: - nrcv = comm.recv( source=r, tag=1 ) - fieldRcv = np.zeros( nrcv, dtype=np.float64 ) - ghostRankRcv = np.zeros( nrcv, dtype=np.int32 ) - localToGlobalMapRcv = np.zeros( nrcv, dtype=np.int64 ) - - comm.Recv( ghostRankRcv, r, 3 ) - comm.Recv( localToGlobalMapRcv, r, 4 ) - - jj = np.where(ghostRankRcv < 0)[0] - fieldRcv[jj] = fullField[localToGlobalMapRcv[jj]] - - comm.Send( fieldRcv, dest=r, tag=100+r ) - - else: - nrcv = nlocalElements - comm.send( nrcv, root, 1 ) - comm.Send( ghostRank, root, 3 ) - comm.Send( localToGlobalMap, root, 4 ) - - comm.Recv( field, source=root, tag=100+rank ) - - return field - - - def filterGhostRank(self, field, **kwargs): - """ - Filter the ghost rank from a GEOS field - - Parameters - ----------- - field : numpy array - Field to filter - - Returns - ------- - field : numpy array - Filtered field - """ - ghostRank = self.getGhostRank(**kwargs) - ind = np.where(ghostRank<0)[0] - - return field[ind] - - - def getWrapper(self, path): - """ - Get the GEOS wrapper - - Parameters - ----------- - path : str - GEOS path - - Returns - -------- - Requested wrapper - """ - if hasattr(self, "solver"): - wrapper = self.solver.get_wrapper(path) - - return wrapper - - - def getGroup(self, path): - """ - Get the GEOS group - - Parameters - ----------- - path : str - GEOS path - - Returns - -------- - Group of the path requested - """ - if hasattr(self, "solver"): - group = self.solver.get_group(path) - - return group diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py deleted file mode 100644 index 34b1d85b..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py +++ /dev/null @@ -1,520 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -from scipy.fftpack import fftfreq, ifft, fft -import numpy as np -import pygeosx - -from .Solver import Solver - -class WaveSolver(Solver): - """ - WaveSolver Object containing methods useful for simulation using wave solvers in GEOS - - Attributes - ----------- - dt : float - Time step for simulation - minTime : float - Min time to consider - maxTime : float - End time to consider - dtSeismo : float - Time step to save pressure for seismic trace - minTimeSim : float - Starting time of simulation - maxTimeSim : float - End Time of simulation - dtWaveField : float - Time step to save fields - sourceType : str - Type of source - sourceFreq : float - Frequency of the source - name : str - Solver name - type : str - Type of solver - Default is None - geosxArgs : GeosxArgs - Object containing GEOSX launching options - geosx : pygeosx instance - - alreadyInitialized : bool - Flag indicating if the initial conditions have been applied yet - firstGeosxInit : bool - Flag for initialize or reinitialize - collectionTargets : list - Output targets for geosx - hdf5Targets : list - HDF5 output targets for geosx - vtkTargets : list - VTK output targets for geosx - """ - def __init__(self, - dt=None, - minTime=0, - maxTime=None, - dtSeismo=None, - dtWaveField=None, - sourceType=None, - sourceFreq=None, - **kwargs): - """ - Parameters - ---------- - dt : float - Time step for simulation - minTime : float - Starting time of simulation - Default is 0 - maxTime : float - End Time of simulation - dtSeismo : float - Time step to save pressure for seismic trace - dtWaveField : float - Time step to save fields - sourceType : str - Type of source - Default is None - sourceFreq : float - Frequency of the source - Default is None - kwargs : keyword args - geosx_argv : list - GEOSX arguments or command line as a splitted line - """ - super().__init__(**kwargs) - - self.name = None - - self.dt = dt - self.minTime = minTime - self.maxTime = maxTime - - self.minTimeSim = minTime - self.maxTimeSim = maxTime - - self.dtSeismo = dtSeismo - - self.sourceType = sourceType - self.sourceFreq = sourceFreq - - self.dtWaveField = dtWaveField - - - def __repr__(self): - string_list = [] - string_list.append("Solver type : " + type(self).__name__ + "\n") - string_list.append("dt : " + str(self.dt) + "\n") - string_list.append("maxTime : " + str(self.maxTime) + "\n") - string_list.append("dtSeismo : " + str(self.dtSeismo) + "\n") - string_list.append("Outputs : " + str(self.hdf5Targets) + "\n" + str(self.vtkTargets) + "\n") - rep = "" - for string in string_list: - rep += string - - return rep - - - def _setTimeVariables(self): - """Initialize the time variables""" - if hasattr(self.xml, "events"): - events = self.xml.events - if self.dt is None: - for event in events["PeriodicEvent"]: - if isinstance(event, dict): - if event["target"] == "/Solvers/" + self.name: - self.dt = float(event['forceDt']) - else: - if event == "target" and events["PeriodicEvent"]["target"] == "/Solvers/" + self.name: - self.dt = float(events["PeriodicEvent"]["forceDt"]) - else: - self.updateTimeVariable("dt") - - if self.maxTime is None: - self.maxTime = float(events["maxTime"]) - self.maxTimeSim = self.maxTime - else: - self.updateTimeVariable("maxTime") - - if self.dtSeismo is None: - if self.type is None: - self._getType() - - solverdict = self.xml.solvers[self.type] - for k, v in solverdict.items(): - if k == "dtSeismoTrace": - self.dtSeismo = float(v) - else: - self.updateTimeVariable("dtSeismo") - - assert None not in (self.dt, self.dtSeismo, self.maxTime) - - - def _setSourceProperties(self): - """Set the frequency and type of source""" - if self.sourceFreq is None: - if self.type is None: - self._getType() - solverdict = self.xml.solvers[self.type] - for k, v in solverdict.items(): - if k == "timeSourceFrequency": - self.sourceFreq = v - - if self.sourceType is None: - if hasattr(self.xml, "events"): - events = self.xml.events - try: - for event in events["PeriodicEvent"]: - if isinstance(event, dict): - if event["target"] == "/Solvers/" + self.name: - self.sourceType = "ricker" + event['rickerOrder'] - else: - if event == "target" and events["PeriodicEvent"]["target"] == "/Solvers/" + self.name: - self.sourceType = "ricker" + events["PeriodicEvent"]["rickerOrder"] - except: - self.sourceType = "ricker2" - - - def initialize(self, rank=0, xml=None): - """ - Initialization or reinitialization of GEOSX - - Parameters - ---------- - rank : int - Process rank - xml : XML - XML object containing parameters for GEOSX initialization. - Only required if not set in the __init__ OR if different from it - """ - super().initialize(rank, xml) - self._setSourceProperties() - - - def updateTimeStep(self, dt): - """ - Update the solver time step - - Parameters - ---------- - dt : double - Time step - """ - self.dt = dt - - - def updateTimeVariable(self, timeVariable): - """ - Overwrite a GEOS time variable - - Parameters - ---------- - timeVariable : str - Name of the time variable to update - """ - if timeVariable == "maxTime": - assert hasattr(self, "maxTime") - self.geosx.get_wrapper("Events/maxTime").value()[0] = self.maxTime - - elif timeVariable == "minTime": - assert hasattr(self, "minTime") - self.geosx.get_wrapper("Events/minTime").value()[0] = self.minTime - - elif timeVariable == "dt": - assert hasattr(self, "dt") - self.geosx.get_wrapper("Events/solverApplications/forceDt").value()[0] = self.dt - - elif timeVariable == "dtSeismo": - assert hasattr(self, "dtSeismo") - self.geosx.get_wrapper("/Solvers/"+self.name+"/dtSeismoTrace").value()[0] = self.dtSeismo - - - def updateSourceFrequency(self, freq): - """ - Overwrite GEOSX source frequency - - Parameters - ---------- - freq : float - Frequency of the source in Hz - """ - self.geosx.get_wrapper("/Solvers/"+self.name+"/timeSourceFrequency").value()[0] = freq - self.sourceFreq = freq - - - - def updateSourceAndReceivers(self, sourcesCoords=[], receiversCoords=[]): - """ - Update sources and receivers positions in GEOS - - Parameters - ---------- - sourcesCoords : list - List of coordinates for the sources - receiversCoords : list - List of coordinates for the receivers - """ - src_pos_geosx = self.solver.get_wrapper( "sourceCoordinates" ).value() - src_pos_geosx.set_access_level( pygeosx.pylvarray.RESIZEABLE ) - - rcv_pos_geosx = self.solver.get_wrapper( "receiverCoordinates" ).value() - rcv_pos_geosx.set_access_level( pygeosx.pylvarray.RESIZEABLE ) - - src_pos_geosx.resize_all(( len( sourcesCoords ), 3 )) - if len( sourcesCoords ) == 0: - src_pos_geosx.to_numpy()[:] = np.zeros(( 0, 3 )) - else: - src_pos_geosx.to_numpy()[:] = sourcesCoords[:] - - - rcv_pos_geosx.resize_all(( len( receiversCoords ), 3 )) - if len( receiversCoords ) == 0: - rcv_pos_geosx.to_numpy()[:] = np.zeros(( 0, 3 )) - else: - rcv_pos_geosx.to_numpy()[:] = receiversCoords[:] - - self.solver.reinit() - - - def evaluateSource(self): - """ - Evaluate source and update on GEOS - Only ricker order {0 - 4} accepted - """ - sourceTypes = ("ricker0", "ricker1", "ricker2", "ricker3", "ricker4") - assert self.sourceType in sourceTypes, f"Only {sourceTypes} are allowed" - - f0 = self.sourceFreq - delay = 1.0 / f0 - alpha = - ( f0 * np.pi )**2 - - nsamples = int( round( ( self.maxTime - self.minTime ) / self.dt )) + 1 - sourceValue = np.zeros(( nsamples, 1 )) - - order = int( self.sourceType[-1] ) - sgn = ( -1 )**( order + 1 ) - - time = self.minTime - for nt in range(nsamples): - - if self.minTime <= - 1.0 / f0: - tmin = - 2.9 / f0 - tmax = 2.9 / f0 - time_d = time - - else: - time_d = time - delay - tmin = 0.0 - tmax = 2.9 / f0 - - if (time > tmin and time < tmax) or ( self.minTime < - 1 / f0 and time == tmin ): - gaussian = np.exp( alpha * time_d**2) - - if order == 0: - sourceValue[nt, 0] = sgn * gaussian - - elif order == 1: - sourceValue[nt, 0] = sgn * ( 2 * alpha * time_d ) * gaussian - - elif order == 2: - sourceValue[nt, 0] = sgn * ( 2 * alpha + 4 * alpha**2 * time_d**2 ) * gaussian - - elif order == 3: - sourceValue[nt, 0] = sgn * ( 12 * alpha**2 * time_d + 8 * alpha**3 * time_d**3 ) * gaussian - - elif order == 4: - sourceValue[nt, 0] = sgn * ( 12 * alpha**2 + 48 * alpha**3 * time_d**2 + 16 * alpha**4 * time_d**4 ) * gaussian - - time += self.dt - - self.updateSourceFrequency(self.sourceFreq) - self.updateSourceValue(sourceValue) - self.sourceValue = sourceValue - - - def updateSourceValue(self, value): - """ - Update the value of the source in GEOS - - Parameters - ---------- - value : array/list - List/array containing the value of the source at each time step - """ - src_value = self.solver.get_wrapper("sourceValue").value() - src_value.set_access_level(pygeosx.pylvarray.RESIZEABLE) - - src_value.resize_all(value.shape) - src_value.to_numpy()[:] = value[:] - - self.maxTimeSim = ( value.shape[0] - 1 ) * self.dt - self.geosx.get_wrapper("Events/minTime").value()[0] = self.minTime - self.sourceValue = value[:] - - - def filterSource(self, fmax): - """ - Filter the source value and give the value to GEOSX. Note that is can also modify the start and end time of simulation to avoid discontinuity. - - Parameters - ----------- - fmax : float/string - Max frequency of the source wanted. The source then have frequencies in the interval [0,fmax+1] - """ - if str(fmax) == "all": - return - - minTime = self.minTime - maxTime = self.maxTime - dt = self.dt - f0 = self.sourceFreq - - sourceValue = self.sourceValue - - pad = int(round(sourceValue.shape[0]/2)) - n = sourceValue.shape[0] + 2 * pad - - tf = fftfreq(n, dt) - y_fft = np.zeros((n,sourceValue.shape[1]), dtype="complex_") - y = np.zeros(y_fft.shape, dtype="complex_") - - for i in range(y_fft.shape[1]): - y_fft[pad:n-pad,i] = sourceValue[:,i] - y_fft[:,i] = fft(y_fft[:,i])# Perform fourier transform - - isup = np.where(tf>=fmax)[0] - imax = np.where(tf[isup]>=fmax+1)[0][0] - i1 = isup[0] - i2 = isup[imax] - - iinf = np.where(tf<=-fmax)[0] - imin = np.where(tf[iinf]<=-fmax-1)[0][-1] - - i3 = iinf[imin] - i4 = iinf[-1] - - for i in range(y_fft.shape[1]): - y_fft[i1:i2,i] = np.cos((isup[0:imax] - i1)/(i2-i1) * np.pi/2)**2 * y_fft[i1:i2,i] - y_fft[i3:i4,i] = np.cos((iinf[imin:-1] - i4)/(i3-i4) * np.pi/2)**2 * y_fft[i3:i4,i] - y_fft[i2:i3,i] = 0 - - for i in range(y.shape[1]): - y[:,i] = ifft(y_fft[:,i])# Perform inverse fourier transform - - it0 = int(round(abs(minTime/dt))) + pad - d = int(round(1/f0/dt)) - - i1 = max(it0 - 4*d, 0) - i2 = int(round(i1 + d/4)) - - i4 = min(n,n - pad + 4*d) - i3 = int(round(i4 - d/4)) - - for i in range(y.shape[1]): - y[i1:i2,i] = np.cos((np.arange(i1,i2) - i2)/(i2-i1) * np.pi/2)**2 * y[i1:i2,i] - y[i3:i4,i] = np.cos((np.arange(i3,i4) - i3)/(i4-i3) * np.pi/2)**2 * y[i3:i4,i] - y[max(i1-d,0):i1,i] = 0.0 - y[i4:min(i4+d,n),i] = 0.0 - - - t = np.arange(minTime-pad*dt, maxTime+pad*dt+dt/2, dt) - - self.updateSourceValue(np.real(y[max(i1-d,0):min(i4+d,n),:])) - self.minTimeSim = t[max(i1-d,0)] - self.maxTimeSim = t[min(i4+d,n-1)] - self.geosx.get_wrapper("Events/minTime").value()[0] = self.minTimeSim - self.sourceValue = np.real(y[max(i1-d,0):min(i4+d,n),:]) - - - def updateVelocityModel(self, vel, velocityName, **kwargs): - """ - Update velocity value in GEOS - - Parameters - ---------- - vel : float/array - Value(s) for velocity field - velocityName : str - Name of the velocity array in GEOS - """ - prefix = self._getPrefixPath(**kwargs) - - velocity = self.solver.get_wrapper(prefix+velocityName).value() - velocity.set_access_level(pygeosx.pylvarray.MODIFIABLE) - - velocity.to_numpy()[vel>0] = vel[vel>0] - - - def updateDensityModel( self, density, densityName, **kwargs ): - """ - Update density values in GEOS - - Parameters - ----------- - density : array - New values for the density - densityName : str - Name of density array in GEOS - """ - prefix = self._getPrefixPath( **kwargs ) - - d = self.solver.get_wrapper( prefix + densityName ).value() - d.set_access_level( pygeosx.pylvarray.MODIFIABLE ) - - d.to_numpy()[density>0] = density[density>0] - - - def outputWaveField(self, time): - """ - Trigger the wavefield output - - Parameters - ---------- - time : float - Current time of simulation - """ - self.collections[0].collect(time, self.dt) - self.hdf5Outputs[0].output(time, self.dt) - - - def getVelocityModel(self, velocityName, filterGhost=False, **kwargs): - """ - Get the velocity values - - Parameters - ----------- - velocityName : str - Name of velocity array in GEOS - filterGhost : bool - Filter the ghost ranks - - Returns - ------- - Numpy Array : Array containing the velocity values - """ - velocity = self.getField(velocityName, **kwargs) - - if filterGhost: - velocity = self.filterGhostRank(velocity, **kwargs) - - return velocity - - - def getWaveField(self): - pass - - def getWaveFieldAtReceivers(self, comm): - pass diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py deleted file mode 100644 index be13e6dc..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py +++ /dev/null @@ -1,30 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -"""Solvers classes""" - -from .Solver import Solver -from .AcousticSolver import AcousticSolver -from .ReservoirSolver import ReservoirSolver -from .ElasticSolver import ElasticSolver -from .WaveSolver import WaveSolver -from .GeomechanicsSolver import GeomechanicsSolver - -from .utils.solverutils import ( - print_group, - print_with_indent, - printGeosx, - printSolver, - printGroup, -) diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py deleted file mode 100644 index 5c0b5183..00000000 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py +++ /dev/null @@ -1,43 +0,0 @@ -# ------------------------------------------------------------------------------------------------------------ -# SPDX-License-Identifier: LGPL-2.1-only -# -# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 TotalEnergies -# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University -# Copyright (c) 2023-2024 Chevron -# Copyright (c) 2019- GEOS/GEOSX Contributors -# Copyright (c) 2019- INRIA project-team Makutu -# All rights reserved -# -# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. -# ------------------------------------------------------------------------------------------------------------ - -def print_group(group, indent=0): - print("{}{}".format(" " * indent, group)) - - indent += 4 - print("{}wrappers:".format(" " * indent)) - - for wrapper in group.wrappers(): - print("{}{}".format(" " * (indent + 4), wrapper)) - print_with_indent(str(wrapper.value(False)), indent + 8) - - print("{}groups:".format(" " * indent)) - - for subgroup in group.groups(): - print_group(subgroup, indent + 4) - - -def print_with_indent(msg, indent): - indent_str = " " * indent - print(indent_str + msg.replace("\n", "\n" + indent_str)) - - -def printGeosx(solver): - print_group(solver.geosx) - -def printSolver(solver): - print_group(solver.solver) - -def printGroup(solver, path): - print_group(solver.solver.get_group(path)) \ No newline at end of file