difficulty implementing diffusion module

[molkjam]

I m trying to implement diffusion within my parcels run using the HYCOM ocean model and am running into a few issues whilst following the examples. I am trying to implement uniform diffusivity using the UniformDiffusionKh module, however the example instructs me to add diffusivity to the fieldset as a constant using the following line:

fieldset.add_constant_field("Kh_zonal", 1, mesh="flat")

yet the documentation suggests there is no 'add_constant_field' class. i have overcome this by using:

fieldset.add_constant("Kh_zonal", 1)
fieldset.add_constant("dres", 0.01)

I hope this is appropriate.

Next, the examples tutorial states that 'The diffusion-step can in this case be computed after or before advection, thus allowing you to chain kernels using the + operator'. However, when I attempt to execute parcels using:

pset.execute(AdvectionRK4 + DiffusionUniformKh,                 # the kernel (which defines how particles move)
                      runtime=timedelta(days=60),    # the total length of the run
                      dt=-timedelta(minutes=1),      # the timestep of the kernel (negative timestep for backwards run)
                      output_file=output_file)

i get the error message:

TypeError: unsupported operand type(s) for +: 'function' and 'function'

Any advice on how I could get this working would be greatly appreciated

[erikvansebille]

Thanks for raising this, @molkjam. For the difference between the example and the documentation, I'm pinging @daanreijnders

For the error that you're getting, that is easy to fix. You will need to convert at least one of the two functions to a Kernel object first, otherwise the + operator is not defined. So change AdvectionRK4 to pset.Kernel(AdvectionRK4), i.e.

pset.execute(pset.Kernel(AdvectionRK4) + DiffusionUniformKh,                 # the kernel (which defines how particles move)
                      runtime=timedelta(days=60),    # the total length of the run
                      dt=-timedelta(minutes=1),      # the timestep of the kernel (negative timestep for backwards run)
                      output_file=output_file)

[molkjam]

Hi Erik, thanks for your quick response. Unfortunately with your fix it chucks out the following error:

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-4-b56bd1fd64fa> in <module>()
      6 
      7 output_file = pset.ParticleFile(name="/{$HOME}/plastics/simulations/test60d.nc", outputdt=timedelta(hours=1)) # the file name and the time step of the outputs
----> 8 pset.execute(pset.Kernel(AdvectionRK4) + DiffusionUniformKh,                 # the kernel (which defines how particles move)
      9                       runtime=timedelta(days=60),    # the total length of the run
     10                       dt=-timedelta(minutes=1),      # the timestep of the kernel (negative timestep for backwards run)

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/kernel.py in __add__(self, kernel)
    455     def __add__(self, kernel):
    456         if not isinstance(kernel, Kernel):
--> 457             kernel = Kernel(self.fieldset, self.ptype, pyfunc=kernel)
    458         return self.merge(kernel)
    459 

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/kernel.py in __init__(self, fieldset, ptype, pyfunc, funcname, funccode, py_ast, funcvars, c_include, delete_cfiles)
    150             kernelgen = KernelGenerator(fieldset, ptype)
    151             kernel_ccode = kernelgen.generate(deepcopy(self.py_ast),
--> 152                                               self.funcvars)
    153             self.field_args = kernelgen.field_args
    154             self.vector_field_args = kernelgen.vector_field_args

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/codegenerator.py in generate(self, py_ast, funcvars)
    428 
    429         # Generate C-code for all nodes in the Python AST
--> 430         self.visit(py_ast)
    431         self.ccode = py_ast.ccode
    432 

/local1/scratch/moja/envs/parcels/lib/python3.6/ast.py in visit(self, node)
    251         method = 'visit_' + node.__class__.__name__
    252         visitor = getattr(self, method, self.generic_visit)
--> 253         return visitor(node)
    254 
    255     def generic_visit(self, node):

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/codegenerator.py in visit_FunctionDef(self, node)
    459         for stmt in node.body:
    460             if not (hasattr(stmt, 'value') and type(stmt.value) is ast.Str):  # ignore docstrings
--> 461                 self.visit(stmt)
    462 
    463         # Create function declaration and argument list

/local1/scratch/moja/envs/parcels/lib/python3.6/ast.py in visit(self, node)
    251         method = 'visit_' + node.__class__.__name__
    252         visitor = getattr(self, method, self.generic_visit)
--> 253         return visitor(node)
    254 
    255     def generic_visit(self, node):

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/codegenerator.py in visit_Assign(self, node)
    571     def visit_Assign(self, node):
    572         self.visit(node.targets[0])
--> 573         self.visit(node.value)
    574         if isinstance(node.value, ast.List):
    575             # Detect in-place initialisation of multi-dimensional arrays

/local1/scratch/moja/envs/parcels/lib/python3.6/ast.py in visit(self, node)
    251         method = 'visit_' + node.__class__.__name__
    252         visitor = getattr(self, method, self.generic_visit)
--> 253         return visitor(node)
    254 
    255     def generic_visit(self, node):

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/codegenerator.py in visit_Call(self, node)
    523         else:
    524             for a in node.args:
--> 525                 self.visit(a)
    526                 if a.ccode == 'parcels_customed_Cfunc_pointer_args':
    527                     pointer_args = True

/local1/scratch/moja/envs/parcels/lib/python3.6/ast.py in visit(self, node)
    251         method = 'visit_' + node.__class__.__name__
    252         visitor = getattr(self, method, self.generic_visit)
--> 253         return visitor(node)
    254 
    255     def generic_visit(self, node):

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/codegenerator.py in visit_BinOp(self, node)
    655         self.visit(node.left)
    656         self.visit(node.op)
--> 657         self.visit(node.right)
    658         if isinstance(node.op, ast.BitXor):
    659             raise RuntimeError("JIT kernels do not support the '^' operator.\n"

/local1/scratch/moja/envs/parcels/lib/python3.6/ast.py in visit(self, node)
    251         method = 'visit_' + node.__class__.__name__
    252         visitor = getattr(self, method, self.generic_visit)
--> 253         return visitor(node)
    254 
    255     def generic_visit(self, node):

/local1/scratch/moja/envs/parcels/lib/python3.6/site-packages/parcels/codegenerator.py in visit_Subscript(self, node)
    643         elif isinstance(node.value, IntrinsicNode):
    644             raise NotImplementedError("Subscript not implemented for object type %s"
--> 645                                       % type(node.value).__name__)
    646         else:
    647             node.ccode = "%s[%s]" % (node.value.ccode, node.slice.ccode)

NotImplementedError: Subscript not implemented for object type ConstNode

May this be because Zh_zonal and Zh_meridional and being added as constants (currently set as 1 for both) not constant fields? if so, can you think of a workaround for this, given the absence of the 'fieldset.add_constant_field' class?

Many thanks,

Molly

[erikvansebille]

I think this is because you made an error in one of your kernels. Did you use a parentheses ([ ]) as if Kh was a Field? You shouldn't add these parentheses if Kh are constants

[molkjam]

I don't think so, I have done the following:

fname = 'uvHYCOM.nc'
filenames = {'U': fname, 'V': fname}
variables = {'U': 'water_u', 'V': 'water_v'}
dimensions = {'lat': 'lat', 'lon': 'lon', 'time': 'time'} # In the GlobCurrent data the dimensions are also called 'lon', 'lat' and 'time'

fieldset = FieldSet.from_netcdf(filenames, variables, dimensions,vmin=-2000)#,allow_time_extrapolation=True)
fieldset.add_constant("Kh_zonal",1)
fieldset.add_constant("Kh_meridional",1)

I'm sure there's something probably glaringly obvious here that i'm getting incorrect

[erikvansebille]

The default Kernel assumes Kh to be Fields, see https://github.com/OceanParcels/parcels/blob/master/parcels/kernels/advectiondiffusion.py#L86-L111

So you should use

fieldset.add_constant_field("Kh_zonal", 1)
fieldset.add_constant_field("Kh_meridional", 1)

[molkjam]

this bring me back to my original issue:

AttributeError: 'FieldSet' object has no attribute 'add_constant_field'

[erikvansebille]

That is then probably because you don't have the most recent version of Parcels. Depending on how old your version of Parcels is, you could try

fieldset.add_field("Kh_zonal", 1, lon=0, lat=0)
fieldset.add_field("Kh_meridional", 1, lon=0, lat=0)

[molkjam]

I am using version 2.2.1. add_field generated the following error: TypeError: add_field() got an unexpected keyword argument 'lon'

[erikvansebille]

Oh sorry, should be

fieldset.add_field(Field("Kh_zonal", 1, lon=0, lat=0))
fieldset.add_field(Field("Kh_meridional", 1, lon=0, lat=0))

[molkjam]

@erikvansebille , thanks your workaround worked great! @daanreijnders I have 2.2.1 and fieldset.add_constant_field("Kh_zonal", 1, mesh="flat") is not being recognised

[daanreijnders]

The tutorial indeed includes features that are not available yet in the latest release (2.2.1 at the moment). This includes the fieldset.add_constant_field() method, which was added in #967. It's a convenience method for the lines Erik just sent. Sorry for the confusion!

[molkjam]

Thank you both for the advice and clearing up any confusion