heat_capacity array does not match dimensions of 2D state variable

[brian-rose]

As of climlab 0.7.3, when you initialize a 2D domain like this:

state = climlab.column_state(num_lev=50, num_lat=60, water_depth=10.)
print(state['Tatm'].shape)
print(state['Tatm'].domain.heat_capacity.shape)

the result is

(60, 50)
(50,)

i.e. the heat_capacity array is not broadcast over the latitude dimension. This is causing some broadcasting problems, and may be connected to the 2D diffusive model issue in brian-rose/ClimateModeling_courseware#5

[brian-rose]

Actually on deeper reflection, the heat capacity is per unit surface area so would always be the same at every point along the latitude axis. So it's not clear whether the array should be broadcast at all.

[duncanwp]

Hi Brian, thanks for this wonderful tool and resource! I've been trying to setup a 1d model with latitudinally varying heat capacity and meridional moist diffusion but can't make it work.

I think it's because the advection diffusion is expecting a constant K but this gets broadcast to the shape of heat_capacity. Is there something obvious I'm missing or is this just not supported at the moment? I'd be happy to try and put together a pull-request if you can point me in the right direction.

[brian-rose]

Hmm, that's a use case that I don't remember ever seeing before. Can you post a code snippet that you represents what you're trying to do, and I can take a look at where it fails.

[duncanwp]

Sure:

model = climlab.EBM_seasonal()
# Add a hemisphericaly asymetric heat-capacity
model.Ts.domain.heat_capacity = model.Ts.domain.heat_capacity * (np.heaviside(-model.domains['Ts'].lat.points[:, np.newaxis], 0.5) + 1.)
d = climlab.dynamics.MeridionalMoistDiffusion(state={'Ts': model.state['Ts']}, **model.param)
model.add_subprocess('Diffusion', d)
model.integrate_years(8)

This throws the following error:

Traceback (most recent call last):
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\IPython\core\interactiveshell.py", line 3331, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-463d4ffa9a3a>", line 3, in <module>
    runfile('C:/Users/duncan/AppData/Roaming/JetBrains/PyCharm2020.2/scratches/scratch_97.py', wdir='C:/Users/duncan/AppData/Roaming/JetBrains/PyCharm2020.2/scratches')
  File "C:\Program Files\JetBrains\PyCharm 2020.1\plugins\python\helpers\pydev\_pydev_bundle\pydev_umd.py", line 197, in runfile
    pydev_imports.execfile(filename, global_vars, local_vars)  # execute the script
  File "C:\Program Files\JetBrains\PyCharm 2020.1\plugins\python\helpers\pydev\_pydev_imps\_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\n", file, 'exec'), glob, loc)
  File "C:/Users/duncan/AppData/Roaming/JetBrains/PyCharm2020.2/scratches/scratch_97.py", line 20, in <module>
    d = climlab.dynamics.MeridionalMoistDiffusion(state={'Ts': model.state['Ts']}, **model.param)
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\meridional_moist_diffusion.py", line 129, in __init__
    super(MeridionalMoistDiffusion, self).__init__(D=D, **kwargs)
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\meridional_heat_diffusion.py", line 74, in __init__
    self.D = D
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\meridional_heat_diffusion.py", line 84, in D
    self._update_diffusivity()
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\meridional_moist_diffusion.py", line 137, in _update_diffusivity
    self.K = self.D / heat_capacity * const.a**2 * (1+f)
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\advection_diffusion.py", line 155, in K
    self._compute_advdiff_matrix()
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\advection_diffusion.py", line 185, in _compute_advdiff_matrix
    self._advdiffTriDiag = adv_diff_numerics.advdiff_tridiag(X=self._Xcenter,
  File "C:\Users\duncan\miniconda3\envs\xarray\lib\site-packages\climlab\dynamics\adv_diff_numerics.py", line 330, in advdiff_tridiag
    tridiag[...,inds_lower[0],inds_lower[1]] = lower_diagonal
ValueError: shape mismatch: value array of shape (90,89) could not be broadcast to indexing result of shape (89,1)

As I say, this might not be the best way to do this. I tried specifying water_depths as an array but that doesn't work either.

[brian-rose]

Ok, I see. I don't think you're missing anything obvious. Seems like the problem arises in the MeridionalHeatDiffusion class, because it accepts a thermal conductivity argument in energetic units, and divides by heat_capacity internally to get a diffusivity in units of m**2/s. But the staggered grid solver means that the diffusivity exists at cell boundaries rather than cell centers, so there's an inconsistency when heat_capacity is not constant.

Fixing this may require some interpolation mechanisms. Or some other deeper dive into the numerics of the diffusion solver.

[duncanwp]

OK, I'll take a look - thanks! Presumably the diffusivity should be zero at the outermost boundaries (at the poles) so that interpolation on the internal bounds should work?

[brian-rose]

There's no requirement that diffusivity be zero at boundaries. The boundary condition is set by the input argument prescribed_flux which sits on the staggered grid (at cell boundaries, including the pole) and defaults to zero. With zero prescribed flux, the staggered grid scheme conserves energy. Actually I think the diffusivity values right at the boundaries are ignored.