# Reg. usage of primary_vars in a python script

Dear Sir,

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?

Thanks
Pavan.

Hi Pavan,

one example is this one:

``````...
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 hope that helps.

Best regards,
Christoph

Dear Christoph,

Thank you for the clarification. I will try this and get back to you.

Pavan.

Dear Christoph,

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.

How to create a source term mesh? Because the overall mesh 1D_ADSD ain’t working.

Thanks
Pavan.

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.

A brief explanation of this is available here.

And you can use the `identifySubdomains` tool of OGS to accomplish the task.

We are working on a set of tutorials. Then the start with OGS will be easier, hopefully.

Just for the sake of completeness: Currently, `bulk_node_ids` is necessary even if you want to apply a source term to the entire simulation domein.

Dear Chleh,

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.

Please guide me to solve this issue.

Thanks
Pavan.

The error message suggests that you pass the “force overwrite flag”. Probably it’s `--force` or `--force-overwrite`, I don’t know by heart.

You can run `identifySubdomains --help` to see a list of all available options.

I hope that helps.

Dear chleh,

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()
``````

Still, I am getting the following error.

Please let me know what the problem is.

Thanks
Pavan

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.