I found that most of the python scripts that try to modify source terms didn’t utilise the primary variables. I would like to give my source term as a function of primary variables using the defined primary_vars command in python.
But, here is the problem. How to call the primary variables. Because, most of the time, they are vectors (i.e. ux, uy or sigma_x, sigma_y, sigma_xy). How to write the script for these, as they might be vectors. Also, how to write the “Jac” term (i.e. Jacobian) for the specific source terms?
...
class BCM_MonitoringTraction(OpenGeoSys.BoundaryCondition):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
pp_actual = primary_vars[0]
ux_actual = primary_vars[1]
uy_actual = primary_vars[2]
uy_target = ExternalDisplacement(x, t)
uy_diff = uy_actual - uy_target
if abs(uy_diff) > 1e-15:
print("WARNING: uy_diff = ", uy_diff)
value = 0.0
jacob = [0.0, 0.0, 0.0]
return (False, value, jacob)
...
It’s not complete in the sense that the Jacobian is all-zero. The Jacobian is the derivative of the returned flux (value in the example) w.r.t. all primary variables.
The first value True/False in return (False, value, jacob) specifies if a BC should be set at that location. You can use this value to switch between Dirichlet and Neumann BCs: A BoundaryCondition subclass can have both a getFlux() and a getDirichletBCValue() method. There are also examples for Dirichlet BCs in the script I mentioned.
I have tried including a python source term in one of the solved 1D example (AdvectionDiffusionSorptionDecay).
try:
import ogs.callbacks as OpenGeoSys
except ModuleNotFoundError:
import OpenGeoSys
import numpy as np
# Choose parametrization
Km = 0.01 # m
# source term for the benchmark
class KreactSourceTerm(OpenGeoSys.SourceTerm):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
Cs = primary_vars[0]
p = primary_vars[1]
value = Cs/(Km+Cs)
Jac = [(Km/((Km+Cs)^2)), 0]
return (value, Jac)
# source term for the benchmark
class PressureSourceTerm(OpenGeoSys.SourceTerm):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
value = 0
Jac = [0, 0]
return (value, Jac)
# instantiate source term object referenced in OpenGeoSys' prj file
react_source_term = KreactSourceTerm()
p_source_term = PressureSourceTerm()
I got the following error regarding the source term mesh.
OGS complains that the nodal property bulk_node_ids is missing in the source term mesh.
This property is necessary to relate the source term mesh to the bulk mesh/simulation domain.
I have tried generating source term mesh. However, the command identifySubdomains ain’t working. I am not able to generate the bulk_node_ids for the entire domain through that command.
The command recognizes the bulk_node_ids in sub meshes and displays “There is already a ‘bulk_element_ids’ property present in the subdomain mesh” and the program stops with a runtime error.
I am not getting how to generate the source term mesh properly.
I have tried including source terms using python script. However, I am facing an issue in forming the jacobian. As per your suggestion earlier, I formed the jacobian as the derivatives of value w.r.t primary_vars
.
# Collection of python boundary condition (BC) classes for OpenGeoSys
try:
import ogs.callbacks as OpenGeoSys
except ModuleNotFoundError:
import OpenGeoSys
import numpy as np
# Choose parametrization
Km = 0.01 # m
Krea = 1e-5
# source term for the benchmark
class KreactSourceTerm(OpenGeoSys.SourceTerm):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
Bf = primary_vars[0] # Concentration
Bs = primary_vars[1] # Concentration
Cu = primary_vars[2] # Concentration
Cnh = primary_vars[3] # Concentration
Cca = primary_vars[4] # Concentration
Ccal = primary_vars[5] # Concentration
p = primary_vars[6] # Pressure
value = Bs*(Krea*Cu)/(Km+Cu) # Source term expression
Jac = [0, (Krea*Cu)/(Km+Cu), ((Bs*Krea*Km)/(((Km+Cu)*(Km+Cu)))), 0, 0, 0, 0] # Jacobian (d(value, primary_var))
return (value, Jac)
# instantiate source term object referenced in OpenGeoSys' prj file
react_source_term = KreactSourceTerm()
As the error message suggests, there is an index out of range in line 23 of your script. Did you check that you have enough primary variables? You can use print statements in the Python script for debugging.
Maybe in the case of staggered coupling only one/a few primary variables are passed to the source term/BC.
