Initialise from Python

Here is a screenshot showing how to set a field from a python function. This is a prescribed field, but the same technique can be used to set the initial condition for a prognostic field. The python option must be selected instead of constant or from file.

Here is the python function again. The function must be a function of X (the location on the mesh) and t (the current time). X[0] gives the x-coordinate for the current position, with X[1] and X[2] giving y and z.

The Python function will be evaluated for each node in the mesh at the start of the simulation to populate the field values. For time-varying prescribed fields, the function will be evaluated again at the beginning of every timestep to update the field value. If it is known that the value of the field does not in fact vary in time, then the re-evaluation of the Python function on each timestep can be inhibited by setting the .../prescribed/do not recalculate option.

 def val(x, t):
  from math import tanh
  alpha=2.0e-4
  rho_0=1025.0
  d_rho=4.0
  rho_c=1024.0
  d_r=50000*((3)**0.5)
  r_h=( (x[0]-100000)**2+(x[1]-100000)**2 )**0.5
  r_v=-1000*x[2]
  rho_h = rho_c - d_rho*tanh(r_h/d_r)
  rho = rho_h*(1.0 - tanh(r_v/d_r)) + rho_0*tanh(r_v/d_r)
  T = (1.0 - (rho/rho_0))/alpha
  return T

For a scalar field, the function must return a single floating point value. Similarly for a vector field, a sequence of values of length equal to the field dimension must be returned. For a tensor field, there are two cases. For an isotropic field specified with .../value/isotropic/python, the function must return a single float which will be used for all the diagonal entries of the tensor at that point. The off-diagonal entries will be set to zero. For the anisotropic case, the function must return a two- dimensional array (a sequence of sequences). It is the user’s responsibility to ensure that the tensor is symmetric in cases where it should be.

Here is a function returning a vector. This function was written to give the positions of an array of detectors.

 def val(t):
   import math
   ret=[]
     for k in range(100,2000,100):
        for j in range(7000,25100,100):
           for i in range(7000,25100,100):
             ret.append([i,j,k])
  return ret

This function returns a two-dimensional solid rotating vector field about the origin.

 def val(X,t):
   return (-X[1],X[0])

A Python function returning a scalar tracer value which is 1.0 in a disk of radius 0.25 about the point (-0.5,0) and zero elsewhere. This example illustrates that it is possible to use a range of Python commands including importing modules.

 def val(X,t):
   from numpy import matrix
   from math import sqrt
   dx= (matrix(X)-matrix((-0.5,0)))
   r=sqrt(dx*dx.T)
   if (r<0.25):
     return 1.0
   else:
     return 0

A Python function returning a tensor which might be a typical three-dimensional diffusivity. The diffusivity ramps up from zero over the first 50 time units:

 def val(X,t):
   t_ramp=50.
   if t<t_ramp:
     scale=t/t_ramp
   else:
     scale=1.
   return [[scale*100, 0, 0], [0, scale*100, 0], [0, 0, scale*1.e-3]]