Skip to content

Commit

Permalink
Merge pull request #63 from deitrr/dev
Browse files Browse the repository at this point in the history 
Dev
  • Loading branch information
@deitrr
deitrr authored Feb 22, 2022
2 parents 378ba84 + 89cdd6d commit 701ef09
Show file tree
Hide file tree
Showing 11 changed files with 205 additions and 105 deletions.
2 changes: 1 addition & 1 deletion ifile/input_template.thr
View file
Original file line number Diff line number Diff line change
Expand Up @@ -230,7 +230,7 @@ Dv_sponge = 0.01
# bottom of diffusive sponge layer (fractional height)
ns_diff_sponge = 0.75

# order of diffusion operator in diff sponge (2 or 4)
# order of diffusion operator in diff sponge (multiple of 2 <= HyDiffOrder)
order_diff_sponge = 2

#-- Radiative transfer (double gray and two streams ) common options -------------------#
Expand Down
20 changes: 12 additions & 8 deletions mjolnir/hamarr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -842,7 +842,6 @@ def regrid(resultsf, simID, ntsi, nts, pgrid_ref='auto', overwrite=False, comp=4
if input.surface:
surf = 1
extrap_low = (grid.Altitudeh[0] - grid.Altitude[1]) / (grid.Altitude[0] - grid.Altitude[1])
Psurf = output.Pressure[:, 1, :] + extrap_low * (output.Pressure[:, 0, :] - output.Pressure[:, 1, :])

else:
surf = 0
Expand Down Expand Up @@ -887,6 +886,8 @@ def regrid(resultsf, simID, ntsi, nts, pgrid_ref='auto', overwrite=False, comp=4
outall = GetOutput(resultsf, simID, t, t, rotation=rotation, theta_z=theta_z, theta_y=theta_y)
output = outall.output
output.load_reshape_all(grid)
if surf:
Psurf = output.Pressure[:, 1, :] + extrap_low * (output.Pressure[:, 0, :] - output.Pressure[:, 1, :])

# interpolate in z to middle of layers using this
interpz = (grid.Altitude - grid.Altitudeh[:-1]) / (grid.Altitudeh[1:] - grid.Altitudeh[:-1])
Expand Down Expand Up @@ -1381,9 +1382,9 @@ def vertical_lat(input, grid, output, rg, sigmaref, z, slice=['default'], save=T
levp = 40
else:
if use_p:
levp = np.arange(np.ceil(np.nanmin(Zonallt[:, prange[0]]) / csp) * csp, np.floor(np.nanmax(Zonallt[:, prange[0]]) / csp) * csp, csp)
levp = np.arange(np.ceil(np.nanmin(Zonallt[:, prange[0]]) / csp) * csp, np.floor(np.nanmax(Zonallt[:, prange[0]]) / csp) * csp + csp, csp)
else:
levp = np.arange(np.ceil(np.nanmin(Zonallt[:, hrange[0]]) / csp) * csp, np.floor(np.nanmax(Zonallt[:, hrange[0]]) / csp) * csp, csp)
levp = np.arange(np.ceil(np.nanmin(Zonallt[:, hrange[0]]) / csp) * csp, np.floor(np.nanmax(Zonallt[:, hrange[0]]) / csp) * csp + csp, csp)

c2 = ax.contour(latp * 180 / np.pi, ycoord, zvals, levels=levp, colors=cover_color, linewidths=1)
ax.clabel(c2, inline=True, fontsize=6, fmt=clabel_format, use_clabeltext=True)
Expand Down Expand Up @@ -1688,7 +1689,8 @@ def vertical_lon(input, grid, output, rg, sigmaref, z, slice=['default'], save=T

def horizontal_lev(input, grid, output, rg, Plev, z, save=True, axis=False,
wind_vectors=False, use_p=True, clevs=[40], cbar = True,
wind_key_loc=[0.85,-0.1], wind_key_max='default',cmap_center=False):
wind_key_loc=[0.85,-0.1], wind_key_max='default',cmap_center=False,
clabel_format='%#.3g',wind_color='0.8'):
# Set the latitude-longitude grid.
loni, lati = np.meshgrid(rg.Longitude[:], rg.Latitude[:])

Expand Down Expand Up @@ -1780,7 +1782,7 @@ def horizontal_lev(input, grid, output, rg, Plev, z, save=True, axis=False,
if len(clevs) == 1:
clevels = np.int(clevs[0])
elif len(clevs) == 3:
clevels = np.linspace(np.int(clevs[0]), np.int(clevs[1]), np.int(clevs[2]))
clevels = np.linspace(clevs[0], clevs[1], np.int(clevs[2]))
else:
raise IOError("clevs not valid!")
# clevels = np.linspace(900,1470,58)
Expand Down Expand Up @@ -1837,16 +1839,16 @@ def horizontal_lev(input, grid, output, rg, Plev, z, save=True, axis=False,
umax = wind_key_max
del Uiii, Viii
if z['llswap']:
q = ax.quiver(latq, lonq, V, U, color='0.8',scale=10*umax)
q = ax.quiver(latq, lonq, V, U, color=wind_color,scale=10*umax)
else:
q = ax.quiver(lonq, latq, U, V, color='0.8',scale=10*umax)
q = ax.quiver(lonq, latq, U, V, color=wind_color,scale=10*umax)
if wind_key_loc:
ax.quiverkey(q, X=wind_key_loc[0], Y=wind_key_loc[1], U=umax, label='%#.2f m/s' % umax, labelpos='E')

divider = make_axes_locatable(ax)
if cbar == True:
cax = divider.append_axes('right', size='5%', pad=0.05)
kwargs = {'format': '%#d'}
kwargs = {'format': clabel_format}
clb = fig.colorbar(C, cax=cax, orientation='vertical', **kwargs)
clb.set_label(z['label'])

Expand Down Expand Up @@ -1961,6 +1963,8 @@ def streamf_moc_plot(input, grid, output, rg, sigmaref, save=True, axis=False,
else:
raise IOError("'axis = {}' but {} is not an axes.SubplotBase instance".format(axis, axis))

print(np.nanmin(sf[:, prange[0]].T/(10**cscale)),np.nanmax(sf[:, prange[0]].T/(10**cscale)))

for cc in C.collections:
cc.set_edgecolor("face") # fixes a stupid bug in matplotlib 2.0

Expand Down
24 changes: 22 additions & 2 deletions mjolnir/mjolnir_plot_helper.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,15 +60,15 @@ def make_plot(args, save=True, axis=None):
'w0prof', 'g0prof', 'spectrum',
'phase','all','eddyKE','eddyMomMerid','eddyTempMerid','eddyTempVar',
'Etotlev','AngMomlev', 'Entropylev','Kdiffprof', 'RiB','BLheight',
'RTbalance','Riprof','RTbalanceTS','tradprof']
'RTbalance','Riprof','RTbalanceTS','tradprof','taulwprof','Fsens']

rg_needed = ['Tver', 'Tlonver', 'uver', 'ulonver', 'vver', 'wver', 'wlonver', 'Tulev', 'PTver', 'PTlonver', 'ulev', 'PVver', 'PVlev',
'RVlev', 'stream', 'tracer', 'Tsurf', 'insol', 'massf', 'pause_rg',
'mustar', 'qheat',
'TSfuptot', 'TSfdowntot', 'TSfnet', 'TSqheat',
'DGfuptot', 'DGfdowntot', 'DGfnet', 'DGqheat',
'all','eddyKE','eddyMomMerid','eddyTempMerid','eddyTempVar',
'Etotlev', 'AngMomlev', 'Entropylev','RiB','BLheight',] # these types need regrid
'Etotlev', 'AngMomlev', 'Entropylev','RiB','BLheight','Fsens'] # these types need regrid

openrg = 0

Expand Down Expand Up @@ -410,6 +410,16 @@ def make_plot(args, save=True, axis=None):
pfile = call_plot('Tsurf',ham.horizontal_lev,input, grid, output, rg, PR_LV, z, wind_vectors=True, use_p=use_p, clevs=args.clevels, save=save, axis=axis)
plots_created.append(pfile)

if ('Fsens' in pview or 'all' in pview) and (input.BL):
if not hasattr(rg, "F_sens"):
raise ValueError("'F_sens' not available in regrid file")
rg.load(['F_sens'])
PR_LV = np.max(output.Pressure) # not important here
z = {'value': rg.F_sens, 'label': r'Surface heat flux (W m$^{-2}$)', 'name': 'Fsens',
'cmap': 'magma', 'lat': rg.Latitude, 'lon': rg.Longitude, 'mt': maketable, 'llswap': args.latlonswap}
pfile = call_plot('Fsens',ham.horizontal_lev,input, grid, output, rg, PR_LV, z, wind_vectors=True, use_p=use_p, clevs=args.clevels, save=save, axis=axis)
plots_created.append(pfile)

if ('DGfuptot' in pview or 'all' in pview) and input.RT:
# Averaged temperature and wind field (longitude vs latitude)
# PR_LV - Pressure level (Pa)
Expand Down Expand Up @@ -802,6 +812,16 @@ def make_plot(args, save=True, axis=None):
pfile = call_plot('Ri',ham.profile,input, grid, output, z, stride=20, save=save, axis=axis, use_p=False)
plots_created.append(pfile)

if ('taulwprof' in pview or 'all' in pview) and input.RT:
output.load_reshape(grid,['tau'])
# dz = (grid.Altitudeh[1:]-grid.Altitudeh[:-1])
# tau_lw = np.trapz(output.tau_lw[:,::-1,:],x=grid.Altitudeh[None,::-1,None],axis=1)
tau_lw = np.cumsum(output.tau_lw[:,::-1,:],axis=1)
z = {'value': tau_lw[:,::-1,:], 'label': r'$\tau_{LW}$', 'name': 'tau_lw'}
pfile = call_plot('tau_lw',ham.profile,input, grid, output, z, stride=20, save=save, axis=axis)
plots_created.append(pfile)


# --- Global diagnostics -----------------------------------
if 'cons' in pview: # RD: needs some work!
ham.conservation(input, grid, output)
Expand Down
20 changes: 13 additions & 7 deletions src/ESP/esp_initial.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,7 @@ __host__ ESP::ESP(int * point_local_,
bool output_mean,
bool out_interm_momentum,
bool output_diffusion,
bool DiffSponge,
init_PT_profile_types init_PT_profile_,
double Tint_,
double kappa_lw_,
Expand Down Expand Up @@ -205,11 +206,14 @@ __host__ ESP::ESP(int * point_local_,
Csurf = Csurf_config;
//
// Allocate Data
alloc_data(globdiag, output_mean, out_interm_momentum, output_diffusion);
alloc_data(globdiag, output_mean, out_interm_momentum, output_diffusion, DiffSponge);
}

__host__ void
ESP::alloc_data(bool globdiag, bool output_mean, bool out_interm_momentum, bool output_diffusion) {
__host__ void ESP::alloc_data(bool globdiag,
bool output_mean,
bool out_interm_momentum,
bool output_diffusion,
bool DiffSponge) {

//
// Description:
Expand Down Expand Up @@ -398,6 +402,8 @@ ESP::alloc_data(bool globdiag, bool output_mean, bool out_interm_momentum, bool
cudaMalloc((void **)&diffw_d, nv * point_num * sizeof(double));
cudaMalloc((void **)&diffrh_d, nv * point_num * sizeof(double));
cudaMalloc((void **)&diff_d, 6 * nv * point_num * sizeof(double));
if (DiffSponge)
cudaMalloc((void **)&diff_sponge_d, 6 * nv * point_num * sizeof(double));
cudaMalloc((void **)&divg_Mh_d, 3 * nv * point_num * sizeof(double));

cudaMalloc((void **)&Kdh2_d, nv * sizeof(double));
Expand Down Expand Up @@ -1271,11 +1277,11 @@ __host__ bool ESP::initial_values(const std::string &initial_conditions_filename
else {
ksponge = 0;
}
if (order_diff_sponge == 2) {
Kdh2_h[lev] = ksponge * pow(dbar, 2.) / timestep_dyn;
if (order_diff_sponge == sim.HyDiffOrder) {
Kdh4_h[lev] += ksponge * pow(dbar, 1.0 * order_diff_sponge) / timestep_dyn;
}
else if (order_diff_sponge == 4) {
Kdh4_h[lev] += ksponge * pow(dbar, 4.) / timestep_dyn;
else {
Kdh2_h[lev] = ksponge * pow(dbar, 1.0 * order_diff_sponge) / timestep_dyn;
}
}
}
Expand Down
1 change: 0 additions & 1 deletion src/ESP/grid.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2056,7 +2056,6 @@ void Icogrid::set_altitudes_softplus(double *Altitude,
for (int lev = 0; lev < nv; lev++) { //centers of layers
Altitude[lev] = (Altitudeh[lev] + Altitudeh[lev + 1]) / 2.0;
}
printf("stop here");
}


Expand Down
27 changes: 14 additions & 13 deletions src/ESP/thor_driver.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -429,18 +429,15 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
cudaMemset(diff_d, 0, sizeof(double) * 6 * point_num * nv);

cudaDeviceSynchronize();
bool firststep;
int stepnum = 0;
for (int ihyp = 0; ihyp < sim.HyDiffOrder / 2 - 1; ihyp++) {
if (ihyp == 0)
firststep = 1;
else
firststep = 0;
//Updates: diffmh_d, diffw_d, diffrh_d, diffpr_d, diff_d
Diffusion_Op<LN, LN><<<NBD, NT>>>(diffmh_d,
diffw_d,
diffrh_d,
diffpr_d,
diff_d,
diff_sponge_d,
Mhk_d,
Rhok_d,
temperature_d,
Expand All @@ -458,8 +455,8 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
Cp_d,
maps_d,
nl_region,
firststep,
0,
stepnum,
false,
sim.DeepModel,
sim.DiffSponge,
order_diff_sponge,
Expand All @@ -474,6 +471,7 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
diffrh_d,
diffpr_d,
diff_d,
diff_sponge_d,
Mhk_d,
Rhok_d,
temperature_d,
Expand All @@ -492,8 +490,8 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
Cp_d,
point_local_d,
point_num,
firststep,
0,
stepnum,
false,
sim.DeepModel,
sim.DiffSponge,
order_diff_sponge,
Expand All @@ -502,6 +500,7 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
energy_equation,
sim.HyDiffOrder);
cudaDeviceSynchronize();
stepnum += 1;
}

BENCH_POINT_I_S(current_step, rk, "Diffusion_Op1", (), ("diff_d"))
Expand All @@ -516,6 +515,7 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
diffrh_d,
diffpr_d,
diff_d,
diff_sponge_d,
Mhk_d,
Rhok_d,
temperature_d,
Expand All @@ -533,8 +533,8 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
Cp_d,
maps_d,
nl_region,
0,
1,
stepnum,
true,
sim.DeepModel,
sim.DiffSponge,
order_diff_sponge,
Expand All @@ -548,6 +548,7 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
diffrh_d,
diffpr_d,
diff_d,
diff_sponge_d,
Mhk_d,
Rhok_d,
temperature_d,
Expand All @@ -566,8 +567,8 @@ __host__ void ESP::Thor(const SimulationSetup &sim, kernel_diagnostics &diag) {
Cp_d,
point_local_d,
point_num,
0,
1,
stepnum,
true,
sim.DeepModel,
sim.DiffSponge,
order_diff_sponge,
Expand Down
9 changes: 4 additions & 5 deletions src/esp.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -765,11 +765,9 @@ int main(int argc, char** argv) {
}

if (sim.DiffSponge) {
if (order_diff_sponge == 2 || order_diff_sponge == 4) {
config_OK &= true;
}
else {
log::printf("order_diff_sponge config option can only be 2 or 4\n");
if (order_diff_sponge > sim.HyDiffOrder || order_diff_sponge % 2) {
log::printf(
"order_diff_sponge config option must be <= HyDiffOrder and a multiple of 2\n");
config_OK &= false;
}
}
Expand Down Expand Up @@ -1083,6 +1081,7 @@ int main(int argc, char** argv) {
sim.output_mean,
sim.out_interm_momentum,
sim.output_diffusion,
sim.DiffSponge,
init_PT_profile,
Tint,
kappa_lw,
Expand Down
Loading

0 comments on commit 701ef09

Please sign in to comment.