Forcing the differential rotation in spherical shell IVP

[Suprabha Mukhopadhyay]

Dear all,

I was trying to force the differential rotation profile using a forcing term in the Navier-Stokes equation for IVP in a spherical shell. Following is the code snippet with which I tried to do so. 

###

basis = de.ShellBasis(coords, shape=(Nphi, Ntheta, Nr), radii=(Ri, Ro), dealias=dealias, dtype=dtype)

ep = dist.VectorField(coords, bases=basis, name='ep')
ep['g'][0] = 1

...

u_phi_mean = de.Average(u@ep, 'phi')
u_phi_target= dist.Field(bases=basis)
u_phi_target['g'] = r*np.sin(theta)* diff_rot(r,theta)/Ekman
N_rot_cool = 50
tau_cool = N_rot_cool*Ekman
forcing = -(u_phi_mean - u_phi_target) / tau_cool

...

problem.add_equation("dt(u) + grad(p) - viscous_terms - Rayleigh/Prandtl*S*g + τ_L/Ekman + lift(τ_u1, -1) + lift(τ_u2, -2) = -(dot(u,e)) + cross(u, f) + forcing*ep") 

###

However, I ran into an error, probably because of the bases:

###

ValueError: operands could not be broadcast together with shapes (0,16,128) (8,16,128) (8,16,128)

###

It would be very helpful if someone could help me fix the error caused by the forcing term. Many thanks in advance.

Cheers,

Sup