Set permeability with uniform gradient distribution

tjlyh777 · November 19, 2022, 2:21pm

Hi all,
I want to realize a uniform gradient initial permeability distribution from fractured area to intact rock in EGS in OGS6. My idea is to set enough different volume physical groups around the injection and production well and set different uniform permeabilities respectively to realize an approximate uniform gradient distribution. However, it seems to cost too much effort. Is there any other easier way to realize it, like a gradient function? Can somebody give me some clues about this problem. I will be very glad.

Regards,
Yuhao Liu

joergbuchwald · November 21, 2022, 2:35pm

Hi!
there are essentially to ways:

I would assign one material domain and use the Function-type parameter: like it is done here: Tests/Data/ThermoHydroMechanics/BGRaCreepAndInitialStressAtIP_AREHS/arehs-salt-THM01_0.prj · master · ogs / ogs · GitLab
or alternatively: Create a cell array in your bulk VTU mesh like it is done here (under 2.): joss_ogs6py_VTUinterface/demo/paper_ogs6py_vtuio.md at master · joergbuchwald/joss_ogs6py_VTUinterface · GitHub
and read it in as a MeshElement property

tjlyh777 · November 21, 2022, 2:53pm

Hi,joerg
It is very kind of you to give these two ways. I will refer to the links provided.

Best,
Yuhao

tjlyh777 · November 26, 2022, 3:40pm

Hi @joergbuchwald ,
I have tried the first way you provided but met some problems, My model is three-dimensional so that the permeability is a 3*3 tensor. If it is set as a Constant-type parameter the command lines are like:

<parameter>
    <name>k02</name>
    <type>Constant</type>
    <values>1e-15 0 0 0 1e-15 0 0 0 1e-15</values>
</parameter>

If I want to set Function-type parameter I think the command line will be:

<parameter>
            <name>k02</name>
            <type>Function</type>
            <expression>exp(-(((x+100)^2+(y-200)^2)^0.5)*0.06)/1e13+1.16e-16</expression>
            <expression>0</expression>
            <expression>0</expression>
            <expression>0</expression>
            <expression>exp(-(((x+100)^2+(y-200)^2)^0.5)*0.06)/1e13+1.16e-16</expression>
            <expression>0</expression>
            <expression>0</expression>
            <expression>0</expression>
            <expression>exp(-(((x+100)^2+(y-200)^2)^0.5)*0.06)/1e13+1.16e-16</expression>
        </parameter>

However, the result reports an error with a parallel computing version:

[4] info: Initialize processes.
[9] info: Initialize processes.
[0] info: Initialize processes.
[2] info: Initialize processes.
[3] info: Initialize processes.
[5] info: Initialize processes.
[6] info: Initialize processes.
[1] info: Initialize processes.
[8] info: Initialize processes.
[7] info: Initialize processes.
[4] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[1] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[3] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[0] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[2] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[5] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[9] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[7] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[6] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[8] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/FunctionParameter.h:119 operator()()
[7] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[6] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[1] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[2] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[9] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[3] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[5] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[4] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[0] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[8] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[2] error: FunctionParameter: The spatial position has to be set by coordinates.
[1] error: FunctionParameter: The spatial position has to be set by coordinates.
[1] info: OGS terminated on 2022-11-26 22:07:38+0800.
[1] error: OGS terminated with error.
[2] info: OGS terminated on 2022-11-26 22:07:38+0800.
[2] error: OGS terminated with error.
[3] error: FunctionParameter: The spatial position has to be set by coordinates.
[5] error: FunctionParameter: The spatial position has to be set by coordinates.
[3] info: OGS terminated on 2022-11-26 22:07:38+0800.
[3] error: OGS terminated with error.
[5] info: OGS terminated on 2022-11-26 22:07:38+0800.
[5] error: OGS terminated with error.
[7] error: FunctionParameter: The spatial position has to be set by coordinates.
[9] error: FunctionParameter: The spatial position has to be set by coordinates.
[7] info: OGS terminated on 2022-11-26 22:07:38+0800.
[7] error: OGS terminated with error.
[9] info: OGS terminated on 2022-11-26 22:07:38+0800.
[9] error: OGS terminated with error.
[6] error: FunctionParameter: The spatial position has to be set by coordinates.
[6] info: OGS terminated on 2022-11-26 22:07:38+0800.
[6] error: OGS terminated with error.
[4] error: FunctionParameter: The spatial position has to be set by coordinates.
[0] error: FunctionParameter: The spatial position has to be set by coordinates.
[0] info: OGS terminated on 2022-11-26 22:07:38+0800.
[4] info: OGS terminated on 2022-11-26 22:07:38+0800.
[4] error: OGS terminated with error.
[0] error: OGS terminated with error.
[8] error: FunctionParameter: The spatial position has to be set by coordinates.
[8] info: OGS terminated on 2022-11-26 22:07:38+0800.
[8] error: OGS terminated with error.

It seems that the function of permeability has some problems in expression but I don’t know what it is…Could you give me some advice? Thank you!
Best,
Yuhao

joergbuchwald · November 29, 2022, 11:01am

It seems to be not yet implemented for the process you used. (It know, it should work at least for THM, TR and TRM I guess) You could open an issue on the gitlab.
But you could use the MeshElement parameter instead. This should still work, especially if your not doing parameter sweeps.

tjlyh777 · November 29, 2022, 11:29am

Hi joerg,
a little upset to hear that this has not been implemented to the process. Actually, I use the TH2M module and I will open an issue in the gitlab.
For the MeshElement parameter, I have read the second link you provided and download the code from github. I have not used the ogs6py as a python interface before so I want to know how to get ogs6py? Do I have to recompile the OGS source code again and change some options using ccmake? Or are there any relevant documentations on the website for me to get some help? I hope you can give me some advice. Thank you!
Best,
Yuhao

joergbuchwald · November 29, 2022, 12:06pm

The interesting part in the notebook was just an example to create an Cell array in the VTU file from a python function depending x,y,z using VTUinterface. You can of course use pyVista, MeshIO or Paraview instead if you know it better. There is no need to recompile or set a special option in cmake. Just create a cell array in the vtu file using the tool of your choice and read it in as MeshElement parameter as done here Tests/Data/Elliptic/nonuniform_bc_SteadyStateDiffusion/inhomogeneous_permeability.prj · master · ogs / ogs · GitLab

tjlyh777 · November 29, 2022, 1:05pm

Thanks for your reply. I am sorry to forget sharing the permeability distribution I want to realize . It is like figure 1 drawn by matlab, which is complex. I have divided the bulk mesh into 6 materials as figure 2 shows in Paraview. If I use Paraview I think “Calculator” is a good choice but it cannot set the distribution by MaterialIDs. So are there any possibilty to realize such complex permeability distribution in Paraview? (Maybe some division and merge operations?) Or I have to use other tools?

Figure 1 Wanted permeability distribution

Figure 2 MaterialIDs

joergbuchwald · November 29, 2022, 3:47pm

I think it is easier to do it in python.
With VTUinterface you need to basically adapt this example:

github.com

joergbuchwald/VTUinterface/blob/master/examples/vtu_create_field.py

'''
This example demonstrates the usage of vtuIO.VTUIO to create point field (arrays)
on the mesh.
'''


import vtuIO	# to read and process (point interpolation) vtu- and pvd-files 


# read file
data=vtuIO.VTUIO("square_1e2_pcs_0_ts_1_t_1.000000.vtu", dim=2)

def get_pressure(x,y,z):
    if x<=0.5:
        return -5.0e3
    else:
        return 5.0e3

data.func_to_field(get_pressure, "p0","p0_field.vtu")

This file has been truncated. show original

With less confusing code as in the notebook before;-)

tjlyh777 · December 2, 2022, 10:44am

Hi @joergbuchwald
I have installed VTUinterface and successfully imported a heterogeneous permeability as an element property. Since my model is 3D so I think the expression should be [kxx, kyy, kzz, kxy, kyz, kxz], 6 components in total, right? Actually I didn’t know permeability is a cell property until now…Thank you very much! But I still have some question. I will appreciate it very much if you could spend some of your precious time to read my following text.
Now I am considering a similar porosity distribution which should match the permeability in my opinion. Porosity is a nodal property so I set “cell=False”(which is default) to get a similar distribution of porosity. In this way I think the MeshNode Parameter is correct for the input. The code is like:

<parameter>
    <name>porosity</name>
    <type>MeshNode</type>
    <field_name>porosity</field_name>
 </parameter>

But an error occurs again:

 critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[2] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[4] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[5] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[7] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[9] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[3] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[1] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[6] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[8] critical: /mnt/e/ogs-petsc-mfront-version/ogs/ParameterLib/MeshNodeParameter.h:51 operator()()
[8] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[7] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[0] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[5] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[9] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[3] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[4] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[6] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[2] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[1] warning: There are still 2 global vectors in use. This might be an indicator of a possible waste of memory.
[5] error: Trying to access a MeshNodeParameter but the node id is not specified.
[3] error: Trying to access a MeshNodeParameter but the node id is not specified.
[5] info: OGS terminated on 2022-12-02 18:12:31+0800.
[5] error: OGS terminated with error.
[0] error: Trying to access a MeshNodeParameter but the node id is not specified.
[3] info: OGS terminated on 2022-12-02 18:12:31+0800.
[3] error: OGS terminated with error.
[0] info: OGS terminated on 2022-12-02 18:12:31+0800.
[0] error: OGS terminated with error.
[9] error: Trying to access a MeshNodeParameter but the node id is not specified.
[9] info: OGS terminated on 2022-12-02 18:12:31+0800.
[8] error: Trying to access a MeshNodeParameter but the node id is not specified.
[9] error: OGS terminated with error.
[8] info: OGS terminated on 2022-12-02 18:12:31+0800.
[8] error: OGS terminated with error.
[7] error: Trying to access a MeshNodeParameter but the node id is not specified.
[7] info: OGS terminated on 2022-12-02 18:12:31+0800.
[7] error: OGS terminated with error.
[4] error: Trying to access a MeshNodeParameter but the node id is not specified.
[1] error: Trying to access a MeshNodeParameter but the node id is not specified.
[1] info: OGS terminated on 2022-12-02 18:12:31+0800.
[1] error: OGS terminated with error.
[4] info: OGS terminated on 2022-12-02 18:12:31+0800.
[4] error: OGS terminated with error.
[6] error: Trying to access a MeshNodeParameter but the node id is not specified.
[6] info: OGS terminated on 2022-12-02 18:12:31+0800.
[6] error: OGS terminated with error.
[2] error: Trying to access a MeshNodeParameter but the node id is not specified.
[2] info: OGS terminated on 2022-12-02 18:12:31+0800.
[2] error: OGS terminated with error.

It seems that node id is lost in my bulk mesh. However, I am sure I added the bulk node id in the bulk mesh and partitioned_node_properties_cfg10.bin file truly exists for the bulk mesh. When I use VTUinterface I noticed that the newly created mesh often has a small memory than the original mesh. Is there any relationship between the lost of memory and the lost of node id? I hope you could give me some advice about the error and the lost of memory of the new mesh.
Best regards,
Yuhao

joergbuchwald · December 2, 2022, 12:20pm

Yes, 6 component for a symetric tensor.
Usually, material properties live on the elements and primary and secondary variables on the nodes.
Therefore, porosity should also be an element property. This is very likely the reason why node id is not specified.

tjlyh777 · December 2, 2022, 1:04pm

Hi
I have set porosity as an element property and the distribution of permeability and porosity has been imported successfully. It is very kind of you to help me out so much. Thank you!
Best,
Yuhao

joergbuchwald · December 2, 2022, 7:51pm

I opened up a MR !4392 which should fix the Function parameter issue in TH2M.