|
| 1 | +"""XBT instrument.""" |
| 2 | + |
| 3 | +from dataclasses import dataclass |
| 4 | +from datetime import timedelta |
| 5 | +from pathlib import Path |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +from parcels import FieldSet, JITParticle, ParticleSet, Variable |
| 9 | + |
| 10 | +from ..spacetime import Spacetime |
| 11 | + |
| 12 | + |
| 13 | +@dataclass |
| 14 | +class XBT: |
| 15 | + """Configuration for a single XBT.""" |
| 16 | + |
| 17 | + spacetime: Spacetime |
| 18 | + min_depth: float |
| 19 | + max_depth: float |
| 20 | + fall_speed: float |
| 21 | + deceleration_coefficient: float |
| 22 | + |
| 23 | + |
| 24 | +_XBTParticle = JITParticle.add_variables( |
| 25 | + [ |
| 26 | + Variable("temperature", dtype=np.float32, initial=np.nan), |
| 27 | + Variable("max_depth", dtype=np.float32), |
| 28 | + Variable("min_depth", dtype=np.float32), |
| 29 | + Variable("fall_speed", dtype=np.float32), |
| 30 | + Variable("deceleration_coefficient", dtype=np.float32), |
| 31 | + ] |
| 32 | +) |
| 33 | + |
| 34 | + |
| 35 | +def _sample_temperature(particle, fieldset, time): |
| 36 | + particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon] |
| 37 | + |
| 38 | + |
| 39 | +def _xbt_cast(particle, fieldset, time): |
| 40 | + particle_ddepth = -particle.fall_speed * particle.dt |
| 41 | + |
| 42 | + # update the fall speed from the quadractic fall-rate equation |
| 43 | + # check https://doi.org/10.5194/os-7-231-2011 |
| 44 | + particle.fall_speed = ( |
| 45 | + particle.fall_speed - 2 * particle.deceleration_coefficient * particle.dt |
| 46 | + ) |
| 47 | + |
| 48 | + # delete particle if depth is exactly max_depth |
| 49 | + if particle.depth == particle.max_depth: |
| 50 | + particle.delete() |
| 51 | + |
| 52 | + # set particle depth to max depth if it's too deep |
| 53 | + if particle.depth + particle_ddepth < particle.max_depth: |
| 54 | + particle_ddepth = particle.max_depth - particle.depth |
| 55 | + |
| 56 | + |
| 57 | +def simulate_xbt( |
| 58 | + fieldset: FieldSet, |
| 59 | + out_path: str | Path, |
| 60 | + xbts: list[XBT], |
| 61 | + outputdt: timedelta, |
| 62 | +) -> None: |
| 63 | + """ |
| 64 | + Use Parcels to simulate a set of XBTs in a fieldset. |
| 65 | +
|
| 66 | + :param fieldset: The fieldset to simulate the XBTs in. |
| 67 | + :param out_path: The path to write the results to. |
| 68 | + :param xbts: A list of XBTs to simulate. |
| 69 | + :param outputdt: Interval which dictates the update frequency of file output during simulation |
| 70 | + :raises ValueError: Whenever provided XBTs, fieldset, are not compatible with this function. |
| 71 | + """ |
| 72 | + DT = 10.0 # dt of XBT simulation integrator |
| 73 | + |
| 74 | + if len(xbts) == 0: |
| 75 | + print( |
| 76 | + "No XBTs provided. Parcels currently crashes when providing an empty particle set, so no XBT simulation will be done and no files will be created." |
| 77 | + ) |
| 78 | + # TODO when Parcels supports it this check can be removed. |
| 79 | + return |
| 80 | + |
| 81 | + fieldset_starttime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[0]) |
| 82 | + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) |
| 83 | + |
| 84 | + # deploy time for all xbts should be later than fieldset start time |
| 85 | + if not all( |
| 86 | + [np.datetime64(xbt.spacetime.time) >= fieldset_starttime for xbt in xbts] |
| 87 | + ): |
| 88 | + raise ValueError("XBT deployed before fieldset starts.") |
| 89 | + |
| 90 | + # depth the xbt will go to. shallowest between xbt max depth and bathymetry. |
| 91 | + max_depths = [ |
| 92 | + max( |
| 93 | + xbt.max_depth, |
| 94 | + fieldset.bathymetry.eval( |
| 95 | + z=0, y=xbt.spacetime.location.lat, x=xbt.spacetime.location.lon, time=0 |
| 96 | + ), |
| 97 | + ) |
| 98 | + for xbt in xbts |
| 99 | + ] |
| 100 | + |
| 101 | + # initial fall speeds |
| 102 | + initial_fall_speeds = [xbt.fall_speed for xbt in xbts] |
| 103 | + |
| 104 | + # XBT depth can not be too shallow, because kernel would break. |
| 105 | + # This shallow is not useful anyway, no need to support. |
| 106 | + for max_depth, fall_speed in zip(max_depths, initial_fall_speeds, strict=False): |
| 107 | + if not max_depth <= -DT * fall_speed: |
| 108 | + raise ValueError( |
| 109 | + f"XBT max_depth or bathymetry shallower than maximum {-DT * fall_speed}" |
| 110 | + ) |
| 111 | + |
| 112 | + # define xbt particles |
| 113 | + xbt_particleset = ParticleSet( |
| 114 | + fieldset=fieldset, |
| 115 | + pclass=_XBTParticle, |
| 116 | + lon=[xbt.spacetime.location.lon for xbt in xbts], |
| 117 | + lat=[xbt.spacetime.location.lat for xbt in xbts], |
| 118 | + depth=[xbt.min_depth for xbt in xbts], |
| 119 | + time=[xbt.spacetime.time for xbt in xbts], |
| 120 | + max_depth=max_depths, |
| 121 | + min_depth=[xbt.min_depth for xbt in xbts], |
| 122 | + fall_speed=[xbt.fall_speed for xbt in xbts], |
| 123 | + ) |
| 124 | + |
| 125 | + # define output file for the simulation |
| 126 | + out_file = xbt_particleset.ParticleFile(name=out_path, outputdt=outputdt) |
| 127 | + |
| 128 | + # execute simulation |
| 129 | + xbt_particleset.execute( |
| 130 | + [_sample_temperature, _xbt_cast], |
| 131 | + endtime=fieldset_endtime, |
| 132 | + dt=DT, |
| 133 | + verbose_progress=False, |
| 134 | + output_file=out_file, |
| 135 | + ) |
| 136 | + |
| 137 | + # there should be no particles left, as they delete themselves when they finish profiling |
| 138 | + if len(xbt_particleset.particledata) != 0: |
| 139 | + raise ValueError( |
| 140 | + "Simulation ended before XBT finished profiling. This most likely means the field time dimension did not match the simulation time span." |
| 141 | + ) |
0 commit comments