joergbuchwald / ogs6py Goto Github PK

Hi,
Currently there is no option in the processes.py submodule to add the calculatesurfaceflux tag to the prj file
https://doxygen.opengeosys.org/d5/d22/ogs_file_param__prj__processes__process__calculatesurfaceflux.html

I suggest to add a function

    def add_surfaceflux(self,**args):
        """
        Add SurfaceFlux

        Parameters
        ----------
        **args : TYPE
            DESCRIPTION.

        Raises
        ------
        KeyError
            DESCRIPTION.

        Returns
        -------
        None.

        """
        self._convertargs(args)

        if "mesh" not in args:
            raise KeyError("No surface mesh for flux analysis assigned")
        if "property_name" not in args:
            raise KeyError("No property name, e.g specific_flux, assigned")
        self.tree['processes']['children']['process']['children'][
                'calculatesurfaceflux'] = self.sflux_vartree
        
        self.sflux_vartree['children']['mesh'] = {
                'tag': 'mesh',
                'text': args['mesh'],
                'attr': {},
                'children': {}
                }
        self.sflux_vartree['children']['property_name'] = {
                'tag': 'property_name',
                'text': args['property_name'],
                'attr': {},
                'children': {}
                }

self.sflux_vartree needs to defined when class Processes is initialised
e.g.

        self.sflux_vartree ={
            'tag': 'calculatesurfaceflux',
            'text': '',
            'attr': {},
            'children': {}
            }

input:
to be implemented

output:

<initial_stress>initial_stress</initial_stress>

It is important to check, whether a simulation converged or not, so this should be returned in a status flag.
However, a non-converged solution should not thrown a RuntimeError and exit the script, since parameter studies may enter ranges of invalid parameter values and still continue.

Great tool!
It would be nice, if in the preamble of example_replace.py would be mentioned that it refers to test.prj generated by example_THM.py.
And in this example porosity is not a medium but a phase property.

Hello,

I'm reviewing your submission for JOSS and having a tough time getting started installing your library. I'd suggest the following:

• Provide a package of your library on PyPi
• Check for an installation of OGS and then raise a descriptive error when the application is missing. For example, immediately after:
  

  ogs6py/ogs6py/ogs.py

  Lines 114 to 117 in 354945b

  	if sys.platform == "win32":
  	ogs_path = os.path.join(ogs_path, "ogs.exe")
  	else:
  	ogs_path = os.path.join(ogs_path, "ogs")

  
  Add a test if the file exists with os.path.isfile, and then point the user to a url containing application installation directions.
• Avoid non-descriptive runtime errors:
  

  ogs6py/ogs6py/ogs.py

  Line 139 in 354945b

  raise RuntimeError

Cheers!
Alex

It would be nice if run_model() would allow for passing additional arguments to the ogs executable call, e.g. for passing an output directory.

Example:

model.run_model(args="-o _out")

It would be nice if there were arguments to run_model() which specify a wrapper and its arguments, e.g. for running with MPI. Example API:

model.run_model(logfile="out.log", wrapper="mpirun", wrapper_args="-np 3")

This would call mpirun -np 3 ogs path/to/project.prj > out.log behind the scenes.

Joerg, thanks for your package. I cannot run it on Windows 10 since the source code is Linux/Mac OS-specific. The error message follows

FileNotFoundError                         Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_5624/1425457894.py in <module>
      1 model.write_input()
----> 2 model.run_model(logfile="out.log")

d:\venv\lib\site-packages\ogs6py\ogs.py in run_model(self, **args)
    475         cmd += f"{self.prjfile} > {self.logfile}"
    476         startt = time.time()
--> 477         returncode = subprocess.run([cmd], shell=True, executable="/bin/bash")
    478         stopt = time.time()
    479         self.exec_time = stopt - startt

input:

model.media.add_property(medium_id="4", phase_type="AqueousLiquid", name="viscosity", type="Function", 
                            expression="exp(-3.719+578.919/(-137.546+temperature))*1e-3")

output:

                        <property>
                            <name>viscosity</name>
                            <type>Function</type>
                        </property>

Hi Joerg,
we want to use the following MPL model with OGS6PY. But it is not considered yet.

model.media.add_property(medium_id="0",                                                                                   
                    name="permeability",
                    type="GasPressureDependentPermeability",
                    initial_permeability="1e-15",
                    a1="0.125",
                    a2="152",
                    pressure_threshold="3.2e6",
                    minimum_permeability="0.0"	,												   
                    maximum_permeability="8.0e-16")	

THis is actually how the model is translated

                <property>
                    <name>permeability</name>
                    <type>GasPressureDependentPermeability</type>
                </property>

and here is how the model should be:

                <property>
                    <name>permeability</name>
                    <type>GasPressureDependentPermeability</type>
                    <initial_permeability>permeability0</initial_permeability>
                    <a1>0.125</a1>
                    <a2>152</a2>
                    <pressure_threshold>3.2e6</pressure_threshold>
                 <minimum_permeability>0.0</minimum_permeability>													   
                    <maximum_permeability>8.0e-16</maximum_permeability>										  
                </property>

Could you please add this model into ogs6py ?
Thanks, Eric

Hi Jörg,
I am trying to use ogs6py to generate a project file for a simulation. In this simulation, I want to apply over 200 thermal source terms that I specify with model.processvars.set_ic(), I just decovered that this function is overwriting previous source terms already generated using the same function, i.e. it is only possible to precribe a single source term in a file. If you precribe multiple source terms, only the latest one specified with model.processvars.set_ic() will be considered

see the following example:

prescribing three source terms with model.processvars.set_ic()

#-----------------------------------------------------------------
# ------------ Emplacement field N° 6 ----------------------------
#-----------------------------------------------------------------

# ------------ field N° 6 Drift N° 1 ----------------------------
model.processvars.add_st(
    process_variable_name="temperature",
    mesh="rangers_fullmodel_0.76mio_consistent_MIDs_set-MaterialIDs-28",
    type="Volumetric",
    parameter="heat_power_HAW"
) 
# ------------ field N° 6 Drift N° 2 ----------------------------
model.processvars.add_st(
    process_variable_name="temperature",
    mesh="rangers_fullmodel_0.76mio_consistent_MIDs_set-MaterialIDs-34",
    type="Volumetric",
    parameter="heat_power_HAW"
) 
# ------------ field N° 6 Drift N° 3 ----------------------------
model.processvars.add_st(
    process_variable_name="temperature",
    mesh="rangers_fullmodel_0.76mio_consistent_MIDs_set-MaterialIDs-40",
    type="Volumetric",
    parameter="heat_power_HAW"
)

generates only the last source term in the project file:

    <process_variables>
        <process_variable>
            <name>temperature</name>
            <components>1</components>
            <order>1</order>
            <initial_condition>T0</initial_condition>
            <boundary_conditions>
                <boundary_condition>
                    <type>Dirichlet</type>
                    <mesh>rangers_fullmodel_0.76mio_consistent_MIDs_boundary_zmax</mesh>
                    <component>0</component>
                    <parameter>T0</parameter>
                </boundary_condition>
                <boundary_condition>
                    <type>Dirichlet</type>
                    <mesh>rangers_fullmodel_0.76mio_consistent_MIDs_boundary_zmin</mesh>
                    <component>0</component>
                    <parameter>T0</parameter>
                </boundary_condition>
            </boundary_conditions>
            <source_terms>
                <source_term>
                    <type>Volumetric</type>
                    <mesh>rangers_fullmodel_0.76mio_consistent_MIDs_set-MaterialIDs-40</mesh>
                    <parameter>heat_power_HAW</parameter>
                </source_term>
            </source_terms>
        </process_variable>
    </process_variables>

I would appreciate if you can take a look at it. The results of your work would be presented next week at the ogs community meeting ;)

staggered_example.zip

• after first coupling iteration the assembly time of the first Newton iteration appears in the coupling iteration
• sometimes the decimal points in dx_relative gets lost, e.g. time step =1, coupling iteration=2, Newton iteration=1, component=0 there appears 2611 instead of 2.611.
• the coupling convergence appears in the order H1, M2, H2, M3, H3, ... i.e. M1 the residuum of the second process in the first coupling iteration is missing

Overall I think this is a pretty good tool for the OGS community, so well done.

I have some minor suggestions:

• You should document how to run the tests in the README, preferably just after the installation instructions. It is always good to know that something has installed and is working as intended. I originally ran tests\test_ogs6py.py from within the tests directory and got a load of missing file errors. It obviously has to be run from the main os6py directory (all the tests passed when I ran it there 😄).
• You might want to mention somewhere that this requires python >= 3.8 (pretty sure that assignment expressions := don't work in versions < 3.8).

Refs openjournals/joss-reviews#3673

input in ogs6py:
to be implemented

output:

    <parameter>
      <mesh>tunnel</mesh>
      <name>Initial_stress</name>
      <type>Function</type>
      <expression>-12.4e6+0.8*4.7e6</expression>
      <expression>-12.7e6+0.8*4.7e6</expression>
      <expression>0</expression>
      <expression>0</expression>
    </parameter>

I have some problems using the method "model.media.add_property" and the method parameter "value", e.g. model.media.add_property(medium_id="0", name="porosity", type="constant", value="42"). It seems that the OGS tag "value" is not written to the OGS project file? What am I doing wrong?

Thanks for the effort to develop such a great tool!
Christian

e.g., for the use in replace_phase_property()

A possible use-case is finding the optimal stabilization parameter for the staggered scheme, currently implemented in the HM process.
The model is attached, in dependence on its coupling scheme parameter should be optimized how many coupling iterations are needed per time step (in average).

The valid range of this parameter is (0.1 ...1), to provoke a non-convergence (to make sure ogs6py continues and provides information about the unsuccessful simulation) a negative value may be entered.

Since there is no a-priori known value for the optimal stabilization parameter in general, this may be a frequent use-case and useful to be included as test/benchmark/example for other users.

use_case_staggered_scheme.zip

Title: OSError When Including XML File in Different Subdirectory

Description:

When using ogs6py for automating project file configurations in a Python script, I encountered an OSError related to the XML include mechanism.

To Reproduce:

1. Use ogs6py in a Python script (parent_dir/notebook.py) to automate the OGS project file configuration.
2. Have a template project file available such as parent_dir/solution_path/template.prj
3. Initialize an OGS instance like so:

   OGS(INPUT_FILE=f"{solution_path}/template.prj", PROJECT_FILE=f"{solution_path}/derived.prj")

4. In template.prj, include an XML file located in the same directory (parent_dir/solution_path/qcurve.xml) using:

   <include file="qcurve.xml" />

5. Run the Python script.

Expected Behavior:

The XML file should be included without any errors.

Actual Behavior:

An OSError occurs:

OSError: Error reading file 'qcurve.xml': failed to load external entity "qcurve.xml"

Additional Context:

• Interestingly, the error does not occur when the XML file is in the parent directory (parent_dir/qcurve.xml).

Is there any specific path resolving behavior for included XML files when using ogs6py?

Staggered schemes output informations about the coupling iterations

info: [time] Solving process #0 took 0.00184468 s in time step #1  coupling iteration #1
info: ------- Checking convergence criterion for coupled solution of process #0 -------
info: Convergence criterion: |dx|=6.4339e+06, |x|=1.2489e+07, |dx|/|x|=5.1517e-01

The convergence information (dx, x, dx/x) should be accessible for all processes (0=H, 1=M in HM simulations), each coupling iteration and each time step.

Note that the coupling convergence is checked after finishing the Newton-iterations for each process.
After the Newton iterations for H converged, the coupling iterations are checked for the convergence in H.
Similarly, after the Newton iterations for M converged, the coupling iterations are checked for the convergence in M.

Note also, that the coupling iterations compare the current result with the previous coupling iteration, i.e. after the first coupling iteration there is no check, since there is no previous result to compare with (this may be important for the context).

First, thanks Joerg for this awesome tool! I was wondering if it would be possible to make optional the writing of the log file when using the model.run_model() method. This is because, for some long running models (with many time steps) the log file quickly becomes huge in size, and for some users this output may not be needed.

Currently, Robin BCs cant be written with ogs6py. I suggest the following change in processvars.py

    def add_bc(self, **args):
        """
        Adds a boundary condition.

        Parameters
        ----------
        process_variable_name : `str`
        type : `str`
        geometrical_set : `str`
        geometry : `str`
        component : `int` or `str`
        parameter : `str`
        bc_object : `str`
        mesh : `str`
        """
        self._convertargs(args)
        if "process_variable_name" not in args:
            raise KeyError("No process variable name specified.")
        if "type" not in args:
            raise KeyError("No type given.")
        if args['process_variable_name'] not in self.tree['process_variables']['children']:
            raise KeyError("You need to set initial condition for that process variable first.")
        if "boundary_conditions" not in self.tree['process_variables']['children'][
                args['process_variable_name']]['children']:
            self.tree['process_variables']['children'][args['process_variable_name']]['children'][
                    'boundary_conditions'] = self.populate_tree('boundary_conditions', children={})
        boundary_conditions = self.tree['process_variables'][
                    'children'][args['process_variable_name']]['children']['boundary_conditions']
        if "geometrical_set" in args:
            if "geometry" not in args:
                raise KeyError("You need to provide a geometry.")
            cpnts = args.get('component','0')
            boundary_conditions['children'][args['geometrical_set']+args['geometry'] +
                    cpnts] = self.populate_tree('boundary_condition', children={})
            boundary_condition = boundary_conditions['children'][
                    args['geometrical_set'] + args['geometry'] + cpnts]
            boundary_condition['children']['type'] = self.populate_tree('type',
                    text=args['type'], children={})
            boundary_condition['children']['geometrical_set'] = self.populate_tree(
                    'geometrical_set', text=args['geometrical_set'], children={})
            boundary_condition['children']['geometry'] = self.populate_tree(
                                'geometry', text=args['geometry'], children={})
            if "parameter" in args:
                if "component" in args:
                    boundary_condition['children']['component'] = self.populate_tree(
                            'component', text=args['component'], children={})
                boundary_condition['children']['parameter'] = self.populate_tree(
                        'parameter', text=args['parameter'], children={})
            elif "bc_object" in args:
                if "component" in args:
                    boundary_condition['children']['component'] = self.populate_tree(
                            'component', text=args['component'], children={})
                boundary_condition['children']['bc_object'] = self.populate_tree(
                        'bc_object', text=args['bc_object'], children={})
            else:
                raise KeyError("Please provide the parameter for Dirichlet \
                                        or Neumann BC/bc_object for Python BC")
        elif "mesh" in args:
            cpnts = args.get('component','0')
            boundary_conditions['children'][args['mesh']+cpnts] = self.populate_tree(
                    'boundary_condition', children={})
            boundary_condition = boundary_conditions['children'][args['mesh']+cpnts]
            boundary_condition['children']['type'] = self.populate_tree('type',
                    text=args['type'], children={})
            boundary_condition['children']['mesh'] = self.populate_tree(
                        'mesh', text=args['mesh'], children={})
            if "parameter" in args:
                if "component" in args:
                    boundary_condition['children']['component'] = self.populate_tree(
                            'component', text=args['component'], children={})
                boundary_condition['children']['parameter'] = self.populate_tree(
                        'parameter', text=args['parameter'], children={})
            elif "bc_object" in args:
                if "component" in args:
                    boundary_condition['children']['component'] = self.populate_tree(
                            'component', text=args['component'], children={})
                boundary_condition['children']['bc_object'] = self.populate_tree(
                        'bc_object', text=args['bc_object'], children={})
            elif "u_0" in args:
                if "alpha" in args:
                    boundary_condition['children']['alpha'] = self.populate_tree(
                            'alpha', text=args['alpha'], children={})
                boundary_condition['children']['u_0'] = self.populate_tree(
                        'u_0', text=args['u_0'], children={})
                
            else:
                raise KeyError("Please provide the parameter for Dirichlet \
                                    or Neumann BC/bc_object for Python BC")
        else:
            raise KeyError("You should provide either a geometrical set \
                                or a mesh to define BC for.")

I either did not understand this line/the ogs6py usage using ogs in a singularity container, or there might be a bug here.
As the code is checking, if singularity is installed, this line could be skipped and L490 could read:

if shutil.which("singularity") is None:

instead. At least, the current implementation does not work for me and I think the reason is given above.

When there is some error in OGS during model.runModel() there is no useful output (just some generic output). It would be very helpful if the OGS stdout/stderr output is written to the console in this case.