Using a kernel to calculate larval growth parameters, running into issues

[TomMBIO]

Question

Question

Good afternoon,

I am working on a parcels simulation that uses interpolated temperature values to calculate the growth of larvae over time. I'm doing this to create a more mechanistic approach to simulating PLDs. I have tried to do this using the code below, but have run into some issues. Can you help me? I have split the code into three "chunks" where I think the issues may lie.

Supporting code/error messages

# Paste your code within this block

###### Code chunk 1: This code creates a funtion that, using development times for a larval species at 3 temperatures from the literature, interpolates a development time for a larva at a specified temperature. **It is encoded in the global environment,** outside of the kernels (which I think causes issues with the JITParticle and its interaction with C -- more on this later)

# Prepare larval physiological data from Sprung (1984).
T_data = np.array([18, 12, 6]) # Experimental temperatures, 18, 12 and 6* C
t_data40 = np.array([457.68, 666.48, 1554.72]) # Developmental times at respective temperatures. As temperature increases, development time decreases due to the energetics of metabolism.

# Sort data for interpolation
sorted_indices = np.argsort(T_data)
T_data = T_data[sorted_indices]
t_data40 = t_data40[sorted_indices]

# Create the interpolation function
interp_func = interp1d(T_data, t_data40, kind='linear')

# Define the development time function
def get_Dt(T):
    """Return interpolated development time (Dt) for a given temperature T."""
    return float(interp_func(T))

# Example:
T=10
get_Dt(T) # Should return a value of 962 hours to full development, at a temperature of 10 degrees C.


###### Code chunk 2: This code designates the various particle variables I intend to use for my calculations

# Creates a "column" within the particle in which to place data
extra_vars = [
    # Variable('age', dtype=np.float32, initial=0, to_write=False),
    Variable('cohort', dtype=np.int32, to_write='once'), # Cohort designation. Tells us the week in which a particle is released
    Variable('on_land', dtype=np.int32, to_write=False), # An important coastal reflection kernel parameter
    Variable('T', dtype=np.float32), # Temperature value at the respective timestep. Interpolated from fieldset "thetao"
    Variable('Dt', dtype=np.float32), # Estimated development time at the given temperature, which is used to calcualte G below
    Variable('G', dtype=np.float32), # Growth parameter. This represents the size increase of a particle per timestep
    Variable('S', dtype=np.float32, initial=70.), # Cumulative size of the particle, starting at 70 microns
]

Particle = JITParticle.add_variables(extra_vars)


###### Code chunk 3: These are my kernels. The first interpolates temperature from the fieldset: it seems to work well. The second uses this temperature to calculate the expected development time (Code chunk 1). This is where I run into my first error. The third uses the expected development time to calculate the growth of the larvae in micrometers for that hour. And finally, the fourth adds this growth to the larva's overall size (so that they grow cumulatively).

def SampleT(particle, fieldset, time):
    """Sample the temperature field"""
    particle.T = fieldset.T[time, particle.depth, particle.lat, particle.lon]

def InterpolateDt(particle, fieldset, time):
    particle.Dt = get_Dt(particle.T)

def CalculateG(particle, fieldset, time):
    particle.G = 230 / particle.Dt #230 is the number of microns required to reach full development, at which point larvae become able to settle.

def GrowS(particle, fieldset, time):
    particle.S += particle.G


#In my pset.execute, I specified the kernels in order of calculation:
kernels = pset.Kernel(Sample_land) + pset.Kernel(Unbeaching) + \
          pset.Kernel(AdvectionRK4) + pset.Kernel(DiffusionUniformKh) + \
          pset.Kernel(SampleT) + pset.Kernel(InterpolateDt) + pset.Kernel(CalculateG) + \
          pset.Kernel(GrowS) + pset.Kernel(CheckOutOfBounds)

I think it's getting there, but whenever I run this code, I get this huge error, the likes of which I do not understand. I believe it has something to do with C getting confused by the use of get_Dt within my kernel.

# Paste your error message within this block
---------------------------------------------------------------------------
CalledProcessError                        Traceback (most recent call last)
File ~\AppData\Roaming\Python\Python312\site-packages\parcels\compilation\codecompiler.py:235, in CCompiler._create_compile_process_(self, cmd, src, log)
    234 try:
--> 235     subprocess.check_call(cmd, stdout=logfile, stderr=logfile)
    236 except OSError:

File ~\anaconda3\envs\parcels\Lib\subprocess.py:413, in check_call(*popenargs, **kwargs)
    412         cmd = popenargs[0]
--> 413     raise CalledProcessError(retcode, cmd)
    414 return 0

CalledProcessError: Command '['gcc', '-Wall', '-fPIC', '-std=gnu11', '-IC:\\Users\\benth\\AppData\\Roaming\\Python\\Python312\\site-packages\\parcels\\include', '-I.', '-g', '-O3', '-DDOUBLE_COORD_VARIABLES', '-m64', '-o', 'C:\\Users\\benth\\AppData\\Local\\Temp\\parcels-tmp\\lib2a058e60beb75c781f32fc12e7081e32_0.dll', 'C:\\Users\\benth\\AppData\\Local\\Temp\\parcels-tmp\\2a058e60beb75c781f32fc12e7081e32_0.c', '-shared', '-lm', '-m64']' returned non-zero exit status 1.

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
Cell In[18], line 29
     23 kernels = pset.Kernel(Sample_land) + pset.Kernel(Unbeaching) + \
     24           pset.Kernel(AdvectionRK4) + pset.Kernel(DiffusionUniformKh) + \
     25           pset.Kernel(SampleT) + pset.Kernel(InterpolateDt) + pset.Kernel(CalculateG) + \
     26           pset.Kernel(GrowS) + pset.Kernel(CheckOutOfBounds)
     28 # Execute simulation
---> 29 pset.execute(
     30     kernels,
     31     runtime=datetime.timedelta(days=60),
     32     dt=datetime.timedelta(minutes=60),
     33     output_file=output_file,
     34 )

File ~\AppData\Roaming\Python\Python312\site-packages\parcels\particleset.py:1053, in ParticleSet.execute(self, pyfunc, pyfunc_inter, endtime, runtime, dt, output_file, verbose_progress, postIterationCallbacks, callbackdt, delete_cfiles)
   1051         self._kernel.remove_lib()
   1052         cppargs = ["-DDOUBLE_COORD_VARIABLES"] if self.particledata.lonlatdepth_dtype else None
-> 1053         self._kernel.compile(
   1054             compiler=GNUCompiler(cppargs=cppargs, incdirs=[os.path.join(get_package_dir(), "include"), "."])
   1055         )
   1056         self._kernel.load_lib()
   1057 if output_file:

File ~\AppData\Roaming\Python\Python312\site-packages\parcels\kernel.py:445, in Kernel.compile(self, compiler)
    442 with open(self.src_file, "w") as f:
    443     f.write(self.ccode)
--> 445 compiler.compile(self.src_file, self.lib_file, self.log_file)
    447 if self.delete_cfiles is False:
    448     logger.info(f"Compiled {self.name} ==> {self.src_file}")

File ~\AppData\Roaming\Python\Python312\site-packages\parcels\compilation\codecompiler.py:311, in GNUCompiler_SS.compile(self, src, obj, log)
    308 lib_pathdir = os.path.dirname(obj)
    309 obj = os.path.join(lib_pathdir, lib_pathfile)
--> 311 super().compile(src, obj, log)

File ~\AppData\Roaming\Python\Python312\site-packages\parcels\compilation\codecompiler.py:275, in CCompiler_SS.compile(self, src, obj, log)
    273 with open(log, "w") as logfile:
    274     logfile.write(f"Compiling: {cc}\n")
--> 275 self._create_compile_process_(cc, src, log)

File ~\AppData\Roaming\Python\Python312\site-packages\parcels\compilation\codecompiler.py:240, in CCompiler._create_compile_process_(self, cmd, src, log)
    238     except subprocess.CalledProcessError:
    239         with open(log) as logfile2:
--> 240             raise RuntimeError(
    241                 f"Error during compilation:\n"
    242                 f"Compilation command: {cmd}\n"
    243                 f"Source/Destination file: {src}\n"
    244                 f"Log file: {logfile.name}\n"
    245                 f"Log output: {logfile2.read()}\n"
    246                 f"\n"
    247                 f"If you are on macOS, it might help to type 'export CC=gcc'"
    248             )
    249 return True

RuntimeError: Error during compilation:
Compilation command: ['gcc', '-Wall', '-fPIC', '-std=gnu11', '-IC:\\Users\\benth\\AppData\\Roaming\\Python\\Python312\\site-packages\\parcels\\include', '-I.', '-g', '-O3', '-DDOUBLE_COORD_VARIABLES', '-m64', '-o', 'C:\\Users\\benth\\AppData\\Local\\Temp\\parcels-tmp\\lib2a058e60beb75c781f32fc12e7081e32_0.dll', 'C:\\Users\\benth\\AppData\\Local\\Temp\\parcels-tmp\\2a058e60beb75c781f32fc12e7081e32_0.c', '-shared', '-lm', '-m64']
Source/Destination file: C:\Users\benth\AppData\Local\Temp\parcels-tmp\2a058e60beb75c781f32fc12e7081e32_0.c
Log file: C:\Users\benth\AppData\Local\Temp\parcels-tmp\2a058e60beb75c781f32fc12e7081e32_0.log
Log output: C:\Users\benth\AppData\Local\Temp\parcels-tmp\2a058e60beb75c781f32fc12e7081e32_0.c:1:0: warning: -fPIC ignored for target (all code is position independent)
 #include "parcels.h"
 ^
C:\Users\benth\AppData\Local\Temp\parcels-tmp\2a058e60beb75c781f32fc12e7081e32_0.c: In function 'Sample_landUnbeachingAdvectionRK4DiffusionUniformKhSampleTInterpolateDtCalculateGGrowSCheckOutOfBounds':
C:\Users\benth\AppData\Local\Temp\parcels-tmp\2a058e60beb75c781f32fc12e7081e32_0.c:198:25: warning: implicit declaration of function 'get_Dt' [-Wimplicit-function-declaration]
   particles->Dt[pnum] = get_Dt(particles->T[pnum]);
                         ^
C:\Users\benth\AppData\Local\Temp\ccbWnPmc.o: In function `Sample_landUnbeachingAdvectionRK4DiffusionUniformKhSampleTInterpolateDtCalculateGGrowSCheckOutOfBounds':
C:/Users/benth/AppData/Local/Temp/parcels-tmp/2a058e60beb75c781f32fc12e7081e32_0.c:198: undefined reference to `get_Dt'
collect2.exe: error: ld returned 1 exit status


If you are on macOS, it might help to type 'export CC=gcc'

Are you able to offer any support with this? I have considered using a non-JIT particle to see if I can get around the C issue. Alternatively, I could try and define the function within the kernel?

Thanks in advance,

Tom

[erikvansebille]

Thanks @TomMBIO, for raising this question. Indeed, your get_DT function is the issue. Kernels in Parcels have a very strict signature - the input parameters (particle, fieldset, time) - which has to be obeyed. So you will need to change your code to something like (not tested)

def get_Dt(particle, fieldset, time):
    # your Kernel code

However, note that in JIT mode, it is not possible to call other functions. Nor can Kernels return a value. The amount of allowed declarations is very limited, see the Kernel part of the Parcels overview tutorial. So you will need to completely change your kernel if you want to run in JIT mode

Note that this limited amount of allowed declarations is actually one of the reasons we're currently working on Parcels v4, but I guess you can't wait another few months until that lands?

[TomMBIO]

Hi Erik, thanks for your timely response.

After writing this I found a work-around by making all the functionality of the development time interpolation function compatible with C:

def CalculateG(particle, fieldset, time): #Where G is the amount a plankton will be expected to grow in an hour, at interpolated temperature (Temp).
    Temp = particle.Temp
    max_growth = 230.0 #Maximum growth before larvae are settlement competent

#Linear interpolation 
    if Temp <= 6.0:
        Dt = 1554.72
    elif Temp >= 18.0:
        Dt = 457.68
    elif Temp < 12.0:
        Dt = 1554.72 + (666.48 - 1554.72) * (Temp - 6.0) / (12.0 - 6.0)
    else:
        Dt = 666.48 + (457.68 - 666.48) * (Temp - 12.0) / (18.0 - 12.0)

    particle.Dt = Dt
    particle.G = max_growth / Dt

This worked well: development time was correctly interpolated, and so was expected growth/hr:
Obs 1 2 3
Temp (*C) 9.405759 9.406226 9.408251
Dt (hours) 1050.5315 1050.4623 1050.1626
G (um) 0.2189368 0.21895123 0.21901372

However, I now have another problem, in that I am struggling to add the value of G (the amount the plankton growns in um per iteration) to the size (S) contained inthe previous iteration:

S (um) 0.2189368 0.21895123 0.21901372 (should be 70 + 0.2189368 0.21895123 0.21901372)

In other words, I want the size (S) of an individual trajectory to be cumulative in much the same way growth is. So far I've written this, which ... didn't work. Any ideas?

def GrowS(particle, fieldset, time):
    if particle.S == 0:
        # First iteration: initialize S to 70 + current G. 70 is the starting size of the larvae in um of the species in question, after fertilisation
        particle.S = 70 + particle.G
    else:
        # Subsequent iterations: accumulate G
        particle.S += particle.G

Also I can't wait that long unfortunately as I need to get this method wrapped up ASAP, but will be interested to see what you add in the next version!

Many thanks,

Tom

[erikvansebille]

Should the first condition in your GrowS kernel not be if particle.S == 70? Because you initialise it at 70, so it will never be 0...

---if particle.S == 0:
+++if particle.S == 70:

[TomMBIO]

I missed this completely. I've changed the initialisation value now and the code now works !

First 10 Growth values: [0.2189368 0.2189396 0.21888779 0.21885504 0.21873033 0.21864954
0.2186085 0.21865638 0.21868007 0.21894722]
First 10 Size values: [70.21893 70.437874 70.65676 70.87562 71.094345 71.312996 71.5316
71.75026 71.96894 72.18789 ]

Cheers,
T