Hi all,
I have some questions about stress initialization of OGS6. I think that in the reality the initial stress is not homogenous in the whole domain. If I have an external geo-mechanical file with stress distribution from an analytical solution or some other computation results from other numerical simulators like Matlab, how can I import it into OGS6 as the initial stress distribution? I know that a stress distribution from another vtu. file using meshadapt command, but how can it be realized from a different kind of file?
I hope someone can provide some suggestions.
Best,
Yuhao
2 Likes
write .vtu file and run OGS from the matlab
@Xue Hi Xue,
Could you please explain it in detail or provide me with a simple example? Thank you very much. 
Best,
Yuhao
Just for repetition, the initial state will be included in the prj-file by lines like
251 <parameter>
252 <name>initial_p</name>
253 <type>MeshNode</type>
254 <field_name>pressure</field_name>
255 <mesh>results_equilibrium_ts_10_t_315360000000.000000</mesh>
256 </parameter>
...
311 <initial_condition>initial_p</initial_condition>
for details about stress I would refer to the code documentation (search in the examples of MeshNode for sigma or stress).
Anyway your questions seems to be how to set up a vtu-file (domain mesh) with the initial values, which obviously are not results from a previous OGS run.
So my recommendation is to add point- or cell-data to this mesh by python code like that
import numpy as np # for operations and data types
from scipy.interpolate import RegularGridInterpolator # for 3D interpolation
import meshio # read/write meshes
inputfile="mwe_input.vtu" # generateStructuredMesh -e hex --lx 10 --ly 10 --lz 10 --nx 20 --ny 20 --nz 20 -o mwe_input.vtu
outputfile="mwe_output.vtu"
nnn = np.array([25, 25, 25]) # divisions of data mesh (to be appended to input mesh) in x-, y-, z-direction
xyz = np.array([0,10, 0,10, 0,10]) # start end (x0,x1, y0,y1, z0,z1)
my_data = np.random.standard_normal(size=(nnn[0], nnn[1], nnn[2])) # generate some data
# data in domain: [X1R,X2R] [Y1R,Y2R] [ZBAS,Z2R]
NXW = nnn[0]
NYW = nnn[1]
NZW = nnn[2]
x = np.linspace(0.0, xyz[1]-xyz[0], NXW) # shifted to x=0 as mesh does
y = np.linspace(0.0, xyz[3]-xyz[2], NYW) # shifted to y=0 as mesh does
z = np.linspace(0.0, xyz[5]-xyz[4], NZW) # shifted to z=0 as mesh does
print("Data grid nxw, nyw, nzw: ", NXW, NYW, NZW)
print("Data grid x1,x2, y1,y2, zb,zr: ", xyz)
# read input mesh
mesh = meshio.read(inputfile)
print(f'Input mesh: {len(mesh.points)} points') # array
print(f'Input mesh: {len(mesh.cells[0].data)} cells') # dictionary
#print(mesh)
# interpolated data to input mesh ([0,LX] [0,LY] [0,LZ])
interpolating_function = RegularGridInterpolator((x, y, z),
my_data[0:NXW, 0:NYW, 0:NZW],
method='nearest',
bounds_error=False,
fill_value=None) # method='nearest', 'linear'
# TODO there may be multiple data on multiple cell blocks
cells_count = len(mesh.cells[0].data) # MeshIO < 5.0 cells_count = len(mesh.cells[0][1])
centerpoints=np.zeros((cells_count, 3)) # evaluate interpolation at center points
for cell_index, cellblock_cell in enumerate(mesh.cells[0].data): # MeshIO < 5.0 for cell_index, cellblock_cell in enumerate(mesh.cells[0][1]):
centerpoints[cell_index] = np.sum(mesh.points[cellblock_cell], axis=0) / len(cellblock_cell)
my_data_interpolation = interpolating_function(centerpoints)
my_data_min = np.min(my_data_interpolation)
my_data_max = np.max(my_data_interpolation)
print("Data appended to input mesh are within range {}...{}".format(my_data_min, my_data_max))
mesh.cell_data['my_data'] = [np.array(my_data_interpolation)] # append data
# write to file
meshio.write(outputfile, mesh)
Without the interpolation stuff the script may be easier, but I do not know what kind of data you have and where they come from.
1 Like
Hi @dominik-kern ,
Thanks for your kind reply. I know the repetition code in prj-file as shown above and what you said
is just what I want to realize and I have csv-file or txt-file generated form Mtalab actually. Thanks for your python code very much. I will learn it carefully.
Best,
Yuhao
Since its customary to prescribe initial data as MeshNodes, you would have to add Point Data instead of the Cell Data in the script attached above.
Another way would be via functions (if such representation is possible), i.e. use lines like that in the prj file
251 <parameter>
252 <name>initial_p</name>
253 <type>Function</type>
254 <expression>9.810000e+03*max(140.000000-z, 0) + 1.000000e+05</expression>
255 </parameter>
256 <parameter>
where the expressions were generated with a script like that (watch for outputs containing prj)
"""
Calculates the initial displacement and pressure of a fully saturated,
homogenous, isotropic elastic cube in a gravitational field
with incompressible grains, and incompressible fluid (alpha=1).
BC-H: pressure on top (atmospheric), no flux elsewhere
BC-M: stress on top (atmospheric), roller elsewhere
@author: dominik kern
"""
# PARAMETERS (average values, constant)
mu = 1.0e-3 # fluid viscosity
k = 1.2e-15 # permeability
E = 7.9e10 # Youngs modulus (bulk)
nu = 0.27 # Poisson ratio (bulk)
a = 1000.0 # characteristic length (to estimate time scale)
rhof = 1000.0 # fluid density
rhos = 2000.0 # solid density
phi = 0.0015 # porosity
g = 9.81 # gravitational acceleration (abs) assumed in negative z-dir.
pa = 1.0e5 # atmospheric pressure
zr_avg = 140.0 # averaged height of top surface (to simplify domain to a cuboid)
z_bottom = -1000 # assuming plane bottom
h_max =zr_avg-z_bottom # maximal depth
rho = phi*rhof + (1-phi)*rhos # bulk density
M = E*(1-nu)/((1+nu)*(1-2*nu)) # P-wave modulus M = K+(4/3)*G
snu = nu/(1-nu) # common factor for stress relation
# TIMESCALE
cv = k*M/mu # [Verruijt: Soil Mechanics] eq. (15.16)
t_ref = (a**2)/cv
print("reference time t_ref = {:2.3e} s".format(t_ref)) # to guide time step
print("--")
# FLUID PRESSURE
# some test values to check for plausibility
h0 = 0.0*h_max # depth of 1. point (given in meters below top)
h1 = 0.5*h_max # depth of 2. point (given in meters below top)
h2 = 1.0*h_max # depth of 3. point (given in meters below top)
p0 = rhof*g*h0 + pa
p1 = rhof*g*h1 + pa
p2 = rhof*g*h2 + pa
print("depth = {:2.3f} m: p = {:2.3e} Pa".format(h0, p0))
print("depth = {:2.3f} m: p = {:2.3e} Pa".format(h1, p1))
print("depth = {:2.3f} m: p = {:2.3e} Pa".format(h2, p2))
# expression for prj-file
print("prj initial_p: {:e}*max({:f}-z, 0) + {:e} (C++ExprTk)".format(rhof*g, zr_avg, pa)) # max() to avoid negative pressures in points above average level
print("--")
# STRESS AND STRAIN
# total stress sigma_zz, some test values to check for plausibility
S0 = -rho*g*h0 - pa
S1 = -rho*g*h1 - pa
S2 = -rho*g*h2 - pa
# effective stress sigma'_zz, some test values to check for plausibility
Se0= S0 + p0
Se1= S1 + p1
Se2= S2 + p2
# output total stresses sigma_xx and sigma_yy
print("depth = {:2.3f} m: sigma_xx = sigma_yy = {:2.3e} Pa, sigma_zz = {:2.3e} Pa".format(h0, snu*Se0-p0, S0))
print("depth = {:2.3f} m: sigma_xx = sigma_yy = {:2.3e} Pa, sigma_zz = {:2.3e} Pa".format(h1, snu*Se1-p1, S1))
print("depth = {:2.3f} m: sigma_xx = sigma_yy = {:2.3e} Pa, sigma_zz = {:2.3e} Pa".format(h2, snu*Se2-p2, S2))
# output effective stresses sigma'_xx and sigma'_yy
print("depth = {:2.3f} m: sigma'_xx = sigma'_yy = {:2.3e} Pa, sigma'_zz = {:2.3e} Pa".format(h0, snu*Se0, Se0))
print("depth = {:2.3f} m: sigma'_xx = sigma'_yy = {:2.3e} Pa, sigma'_zz = {:2.3e} Pa".format(h1, snu*Se1, Se1))
print("depth = {:2.3f} m: sigma'_xx = sigma'_yy = {:2.3e} Pa, sigma'_zz = {:2.3e} Pa".format(h2, snu*Se2, Se2))
# expression for prj-file, only needed when fixed displacement (Dirichlet) is to be enforced approximately by the corresponding stress (Neumann)
print("( prj sigma_xx=sigma_yy: {:2.3e}*({:2.3f}-z) - {:2.3e} )".format( g*(snu*(rhof-rho) - rhof), zr_avg, pa))
print("--")
# DISPLACEMENT (u in z-direction)
uz = lambda z: 0.5*(rho-rhof)*(g/M)*((zr_avg-z)**2 - h_max**2)
print("depth = {:2.3f} m: u_z = {:2.4f} m".format(h0, uz(zr_avg-h0)))
print("depth = {:2.3f} m: u_z = {:2.4f} m".format(h1, uz(zr_avg-h1)))
print("depth = {:2.3f} m: u_z = {:2.4f} m".format(h2, uz(zr_avg-h2)))
print("prj initial_uz: {:e}*(({}-z)^2-{}^2) (C++ExprTk)".format(0.5*(rho-rhof)*g/M, zr_avg, h_max))
By the way, this script is going to be part of the OGS tools.
may i get your email