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 gives the x-coordinate for the current position, with X and X 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-100000)**2+(x-100000)**2 )**0.5 r_v=-1000*x 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
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,X)
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]]