askill
/
claude
mirror of https://github.com/Askill/claude.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Attempted fix of convection 

created pressure background state
This commit is contained in:
Simon Clark 2020-09-30 20:53:03 +01:00
parent 663573f911
commit 54180b050a
2 changed files with 18 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
claude_top_level_library.pyx
View File

@ -53,7 +53,7 @@ def divergence_with_scalar(np.ndarray a,np.ndarray u,np.ndarray v,np.ndarray w,n
for i in range(nlat): for i in range(nlat):
for j in range(nlon): for j in range(nlon):
for k in range(nlevels): for k in range(nlevels):
output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + 1*low_level.scalar_gradient_z(aw,dz,i,j,k) output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + low_level.scalar_gradient_z(aw,dz,i,j,k)
return output return output
# specifically thermal advection, converts temperature field to potential temperature field, then advects (adiabatically), then converts back to temperature # specifically thermal advection, converts temperature field to potential temperature field, then advects (adiabatically), then converts back to temperature
@ -74,7 +74,7 @@ def thermal_advection(np.ndarray temperature_atmos,np.ndarray air_pressure,np.nd
for i in range(nlat): for i in range(nlat):
for j in range(nlon): for j in range(nlon):
for k in range(nlevels): for k in range(nlevels):
output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + 1*low_level.scalar_gradient_z(aw,dz,i,j,k) output[i,j,k] = low_level.scalar_gradient_x(au,dx,nlon,i,j,k) + low_level.scalar_gradient_y(av,dy,nlat,i,j,k) + low_level.scalar_gradient_z(aw,dz,i,j,k)
output = low_level.theta_to_t(output,air_pressure) output = low_level.theta_to_t(output,air_pressure)
@ -129,7 +129,7 @@ def radiation_calculation(np.ndarray temperature_world, np.ndarray temperature_a
return temperature_world, temperature_atmos return temperature_world, temperature_atmos
def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray air_pressure,np.ndarray old_pressure,np.ndarray air_density,np.ndarray coriolis,DTYPE_f gravity,np.ndarray dx,DTYPE_f dy,DTYPE_f dt): def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray w,np.ndarray air_pressure,np.ndarray ref_pressure_profile,np.ndarray air_density,np.ndarray coriolis,DTYPE_f gravity,np.ndarray dx,DTYPE_f dy,np.ndarray dz,DTYPE_f dt):
# introduce temporary arrays to update velocity in the atmosphere # introduce temporary arrays to update velocity in the atmosphere
cdef np.ndarray u_temp = np.zeros_like(u) cdef np.ndarray u_temp = np.zeros_like(u)
cdef np.ndarray v_temp = np.zeros_like(v) cdef np.ndarray v_temp = np.zeros_like(v)
@ -144,22 +144,20 @@ def velocity_calculation(np.ndarray u,np.ndarray v,np.ndarray air_pressure,np.nd
inv_dt_grav = 1/(dt*gravity) inv_dt_grav = 1/(dt*gravity)
# calculate acceleration of atmosphere using primitive equations on beta-plane # calculate acceleration of atmosphere using primitive equations on beta-plane
for i in np.arange(3,nlat-3).tolist(): for i in np.arange(2,nlat-2).tolist():
for j in range(nlon): for j in range(nlon):
for k in range(nlevels): for k in range(nlevels):
inv_density = 1/air_density[i,j,k] inv_density = 1/air_density[i,j,k]
u_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(u,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(u,dy,nlat,i,j,k) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)*inv_density ) u_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(u,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(u,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z(u,dz,i,j,k) + coriolis[i]*v[i,j,k] - low_level.scalar_gradient_x(air_pressure,dx,nlon,i,j,k)*inv_density - 1E-6*u[i,j,k])
v_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(v,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(v,dy,nlat,i,j,k) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(air_pressure,dy,nlat,i,j,k)*inv_density ) v_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(v,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(v,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z(v,dz,i,j,k) - coriolis[i]*u[i,j,k] - low_level.scalar_gradient_y(air_pressure,dy,nlat,i,j,k)*inv_density - 1E-6*v[i,j,k])
w_temp[i,j,k] += -(air_pressure[i,j,k]-old_pressure[i,j,k])*inv_density*inv_dt_grav # w_temp[i,j,k] += dt*( -u[i,j,k]*low_level.scalar_gradient_x(w,dx,nlon,i,j,k) - v[i,j,k]*low_level.scalar_gradient_y(w,dy,nlat,i,j,k) - w[i,j,k]*low_level.scalar_gradient_z(w,dz,i,j,k) - inv_density*low_level.scalar_gradient_z(air_pressure,dz,i,j,k) - gravity - 1E-6*w[i,j,k])
w_temp[i,j,k] += dt*( - inv_density*low_level.scalar_gradient_z_1D((air_pressure[i,j,:]-ref_pressure_profile),dz,k) - gravity - 1E-6*w[i,j,k])
u += u_temp u += u_temp
v += v_temp v += v_temp
w = w_temp w += w_temp
# approximate friction
u *= 0.95
v *= 0.95
# approximate surface friction
u[:,:,0] *= 0.8 u[:,:,0] *= 0.8
v[:,:,0] *= 0.8 v[:,:,0] *= 0.8

16
toy_model.py
View File

@ -22,8 +22,8 @@ axial_tilt = -23.5/2 # tilt of rotational axis w.r.t. solar plane
year = 365*day # length of year (s) year = 365*day # length of year (s)
dt_spinup = 60*137 dt_spinup = 60*137
dt_main = 60*7 dt_main = 60*0.1
spinup_length = 7*day spinup_length = 10*day
### ###
@ -32,7 +32,7 @@ advection_boundary = 3 # how many gridpoints away from poles to apply advectio
smoothing_parameter_t = 0.6 smoothing_parameter_t = 0.6
smoothing_parameter_u = 0.6 smoothing_parameter_u = 0.6
smoothing_parameter_v = 0.6 smoothing_parameter_v = 0.6
smoothing_parameter_w = 0.3 smoothing_parameter_w = 0.2
save = False # save current state to file? save = False # save current state to file?
load = True # load initial state from file? load = True # load initial state from file?
@ -103,6 +103,7 @@ for i in range(nlat):
specific_gas = 287 specific_gas = 287
thermal_diffusivity_roc = 1.5E-6 thermal_diffusivity_roc = 1.5E-6
ref_pressure_profile = density_profile*specific_gas*temp_profile
air_pressure = air_density*specific_gas*temperature_atmos air_pressure = air_density*specific_gas*temperature_atmos
# define planet size and various geometric constants # define planet size and various geometric constants
@ -201,10 +202,10 @@ while True:
if velocity: if velocity:
before_velocity = time.time() before_velocity = time.time()
u,v,w = top_level.velocity_calculation(u,v,air_pressure,old_pressure,air_density,coriolis,gravity,dx,dy,dt) u,v,w = top_level.velocity_calculation(u,v,w,air_pressure,ref_pressure_profile,air_density,coriolis,gravity,dx,dy,dz,dt)
u = top_level.smoothing_3D(u,smoothing_parameter_u) u = top_level.smoothing_3D(u,smoothing_parameter_u)
v = top_level.smoothing_3D(v,smoothing_parameter_v) v = top_level.smoothing_3D(v,smoothing_parameter_v)
w = top_level.smoothing_3D(w,smoothing_parameter_w,0.3) w = top_level.smoothing_3D(w,smoothing_parameter_w,0.1)
u[(advection_boundary,-advection_boundary-1),:,:] *= 0.5 u[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
v[(advection_boundary,-advection_boundary-1),:,:] *= 0.5 v[(advection_boundary,-advection_boundary-1),:,:] *= 0.5
@ -242,7 +243,6 @@ while True:
# time_taken = float(round(time.time() - before_advection,3)) # time_taken = float(round(time.time() - before_advection,3))
# print('Advection: ',str(time_taken),'s') # print('Advection: ',str(time_taken),'s')
# before_plot = time.time() # before_plot = time.time()
if plot: if plot:
@ -268,8 +268,8 @@ while True:
ax[0].axhline(y=0,color='black',alpha=0.3) ax[0].axhline(y=0,color='black',alpha=0.3)
ax[0].set_xlabel('Longitude') ax[0].set_xlabel('Longitude')
ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic',levels=15) # ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(temperature_atmos,axis=1)), cmap='seismic',levels=15)
# ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(w,axis=1)), cmap='seismic',levels=15) ax[1].contourf(heights_plot, lat_z_plot, np.transpose(np.mean(w,axis=1)), cmap='seismic',levels=15)
ax[1].plot(lat_plot,tropopause_height,color='black',linestyle='--',linewidth=3,alpha=0.5) ax[1].plot(lat_plot,tropopause_height,color='black',linestyle='--',linewidth=3,alpha=0.5)
if velocity: if velocity:
ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8) ax[1].contour(heights_plot,lat_z_plot, np.transpose(np.mean(u,axis=1)), colors='white',levels=20,linewidths=1,alpha=0.8)