From b29eb0e3f83ab310b08da5d4bedcc1269c867809 Mon Sep 17 00:00:00 2001 From: Aleks Novikov Date: Wed, 30 Oct 2024 10:11:23 +0100 Subject: [PATCH 1/2] Python interfaces to geos from Makutu repo --- pygeos-tools/pyproject.toml | 4 + .../pygeos_tools/utilities/input/GeosxArgs.py | 202 +++ .../geos/pygeos_tools/utilities/input/Xml.py | 120 ++ .../pygeos_tools/utilities/input/__init__.py | 18 + .../utilities/mesh/InternalMesh.py | 144 +++ .../utilities/mesh/VtkFieldSpecifications.py | 197 +++ .../pygeos_tools/utilities/mesh/VtkMesh.py | 632 ++++++++++ .../pygeos_tools/utilities/mesh/__init__.py | 19 + .../pygeos_tools/utilities/model/__init__ .py | 33 + .../utilities/model/utils/VtkModel.py | 1096 +++++++++++++++++ .../utilities/model/utils/vtkUtils.py | 488 ++++++++ .../utilities/solvers/AcousticSolver.py | 398 ++++++ .../utilities/solvers/ElasticSolver.py | 319 +++++ .../utilities/solvers/GeomechanicsSolver.py | 104 ++ .../utilities/solvers/ReservoirSolver.py | 216 ++++ .../pygeos_tools/utilities/solvers/Solver.py | 697 +++++++++++ .../utilities/solvers/WaveSolver.py | 520 ++++++++ .../utilities/solvers/__init__.py | 30 + .../utilities/solvers/utils/solverutils.py | 43 + 19 files changed, 5280 insertions(+) create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/mesh/__init__.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/model/__init__ .py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/VtkModel.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/vtkUtils.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/AcousticSolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py create mode 100644 pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py diff --git a/pygeos-tools/pyproject.toml b/pygeos-tools/pyproject.toml index 2a6f5dec..24f846fd 100644 --- a/pygeos-tools/pyproject.toml +++ b/pygeos-tools/pyproject.toml @@ -22,6 +22,10 @@ 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 new file mode 100644 index 00000000..a19454e5 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py @@ -0,0 +1,202 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..8c61aaed --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py @@ -0,0 +1,120 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..d4d08eba --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py @@ -0,0 +1,18 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..56e6f213 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py @@ -0,0 +1,144 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..ad80d579 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py @@ -0,0 +1,197 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..1f7a7651 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py @@ -0,0 +1,632 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..5f707d45 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py @@ -0,0 +1,319 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..7e21fb46 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py @@ -0,0 +1,104 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..a6489d7f --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py @@ -0,0 +1,216 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..a5241321 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py @@ -0,0 +1,697 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..b22f5a61 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py @@ -0,0 +1,520 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..2e33aa04 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py @@ -0,0 +1,30 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 new file mode 100644 index 00000000..9a207f21 --- /dev/null +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py @@ -0,0 +1,43 @@ +# ------------------------------------------------------------------------------------------------------------ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 Total, S.A +# 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 From 1a0b4ebaefa7d97294d3bb380c4a66ef91ac54e9 Mon Sep 17 00:00:00 2001 From: Aleks Novikov Date: Fri, 15 Nov 2024 10:49:50 +0100 Subject: [PATCH 2/2] Fix company name in copyright preambles --- pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py | 2 +- pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py | 2 +- pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py | 2 +- .../src/geos/pygeos_tools/utilities/mesh/InternalMesh.py | 2 +- .../geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py | 2 +- pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py | 2 +- pygeos-tools/src/geos/pygeos_tools/utilities/mesh/__init__.py | 2 +- pygeos-tools/src/geos/pygeos_tools/utilities/model/__init__ .py | 2 +- .../src/geos/pygeos_tools/utilities/model/utils/VtkModel.py | 2 +- .../src/geos/pygeos_tools/utilities/model/utils/vtkUtils.py | 2 +- .../src/geos/pygeos_tools/utilities/solvers/AcousticSolver.py | 2 +- .../src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py | 2 +- .../geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py | 2 +- .../src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py | 2 +- pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py | 2 +- .../src/geos/pygeos_tools/utilities/solvers/WaveSolver.py | 2 +- .../src/geos/pygeos_tools/utilities/solvers/__init__.py | 2 +- .../geos/pygeos_tools/utilities/solvers/utils/solverutils.py | 2 +- 18 files changed, 18 insertions(+), 18 deletions(-) diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py b/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py index a19454e5..f6faf41b 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/input/GeosxArgs.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py b/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py index 8c61aaed..161fdab7 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/input/Xml.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py b/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py index d4d08eba..af3b85b4 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/input/__init__.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py index 56e6f213..5019b913 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/InternalMesh.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py index ad80d579..cba7259e 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkFieldSpecifications.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py index 1f7a7651..f89a1e54 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/VtkMesh.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/__init__.py b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/__init__.py index 3c73926d..dd924c2a 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/__init__.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/mesh/__init__.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/model/__init__ .py b/pygeos-tools/src/geos/pygeos_tools/utilities/model/__init__ .py index 11aea12e..13f04511 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/model/__init__ .py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/model/__init__ .py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/VtkModel.py b/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/VtkModel.py index b8f45f1f..b612d82b 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/VtkModel.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/VtkModel.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/vtkUtils.py b/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/vtkUtils.py index ceccd426..47069ff3 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/vtkUtils.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/model/utils/vtkUtils.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/AcousticSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/AcousticSolver.py index b94fe75f..5ef81c1e 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/AcousticSolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/AcousticSolver.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py index 5f707d45..34acbc77 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ElasticSolver.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py index 7e21fb46..ebf16fea 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/GeomechanicsSolver.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py index a6489d7f..48f7a438 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/ReservoirSolver.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py index a5241321..a2afd4d9 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/Solver.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py index b22f5a61..34b1d85b 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/WaveSolver.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 diff --git a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py index 2e33aa04..be13e6dc 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/__init__.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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 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 index 9a207f21..5c0b5183 100644 --- a/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py +++ b/pygeos-tools/src/geos/pygeos_tools/utilities/solvers/utils/solverutils.py @@ -2,7 +2,7 @@ # SPDX-License-Identifier: LGPL-2.1-only # # Copyright (c) 2016-2024 Lawrence Livermore National Security LLC -# Copyright (c) 2018-2024 Total, S.A +# 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