Changes for plotting.

Author: quantheory

.gitignore

@@ -1,3 +1,4 @@
 *.pyf
 *.pyfe
 *.so
+*.eps

build_micro.csh

@@ -46,7 +46,7 @@ set FILES="${DIR}/wv_sat_methods.F90 ${DIR}/micro_mg_utils.F90 ${DIR}/micro_mg2_
 #are separated by spaces - the code still compiles if commas are used instead,
 #but it crashes when you try to load it into python. 
 
-set PUB_FNS='wv_sat_methods_init wv_sat_svp_water wv_sat_svp_ice wv_sat_svp_to_qsat wv_sat_qsat_water wv_sat_qsat_ice micro_mg_init micro_mg_tend calc_precip_frac micro_mg_utils_init ice_deposition_sublimation kk2000_liq_autoconversion ice_autoconversion immersion_freezing contact_freezing snow_self_aggregation accrete_cloud_water_snow secondary_ice_production accrete_rain_snow heterogeneous_rain_freezing accrete_cloud_water_rain self_collection_rain accrete_cloud_ice_snow evaporate_rain sublimate_snow bergeron_process_snow sedimentation avg_diameter rain_fall_speeds'
+set PUB_FNS='wv_sat_methods_init micro_mg_init micro_mg_tend calc_precip_frac micro_mg_utils_init'
 
 #CREATE .pyf FILE WITH INTERFACES TO PYTHON:
 #=================
@@ -99,4 +99,4 @@ sed -i '/type unknown_type/,/type unknown_type/d' ${NM}.pyf
 #Note: -DHAVE_GAMMA_INTRINSICS #defines HAVE_GAMMA_INTRINSICS, which is 
 #true for ifort and prevents an additional awkward dependency.
 
-f2py -c --fcompiler=intelem -DHAVE_GAMMA_INTRINSICS ${NM}.pyf ${FILES}
+f2py -c --fcompiler=gnu95 -DHAVE_GAMMA_INTRINSICS ${NM}.pyf ${FILES}

jacobian2.py

@@ -1,5 +1,7 @@
 #!/usr/bin/env python
 
+import sys
+
 import numpy as np
 import scipy.linalg as la
 import scipy.stats as stats
@@ -14,8 +16,8 @@
 
 from mg2_constants import *
 
-HIST_FILE_NAME = "/g/g14/santos36/Data/MG2_data_collection.cam.h1.0001-01-06-00000.nc"
-CLUSTER_FILE_NAME = "/g/g14/santos36/Data/MG2_data_collection.10_cluster_labels.0001-01-06-00000.nc"
+HIST_FILE_NAME = "/home/santos/Data/MG2_data_collection.cam.h1.0001-01-06-00000.nc"
+CLUSTER_FILE_NAME = "/home/santos/Data/MG2_data_collection.10_cluster_labels.0001-01-06-00000.nc"
 
 hfile = nc4.Dataset(HIST_FILE_NAME, 'r')
 cfile = nc4.Dataset(CLUSTER_FILE_NAME, 'r')
@@ -101,6 +103,12 @@
 frzdep_loc = np.empty((1, frzdep.shape[1]), order='F')
 
 timestep = 1.e-8
+evap_col_steps = 1
+evap_steps = 1
+col_steps = 1
+auto_accr_steps = 1
+auto_steps = 1
+accr_steps = 1
 
 # Indices for each variable in the state array.
 it = 0
@@ -187,6 +195,7 @@ def s_inv_mult(state):
     "iaccr",
     "rnagg",
     "snagg",
+    "cnact",
     "anadj",
 ]
 process_names = {
@@ -229,9 +238,10 @@ def mg2_tendencies(level, t, q, qc, nc, qi, ni, qr, nr, qs, ns):
         nrout2, nsout2, drout2, dsout2, freqs, freqr, nfice, qcrat, \
         nadjtot, ncmeitot, mnudeptot, nnudeptot, npsacwstot, nnuccctot, nnuccttot, \
         nsacwitot, nnuccrtot, mnuccritot, nnuccritot, npracstot, npratot, nprctot, \
-        nprc1tot, npraitot, nprcitot, nsubrtot, nraggtot, nsaggtot, errstring, \
-        prer_evap \
-        = mg.micro_mg_tend(timestep, t[0,level], q[0,level], qc[0,level], qi[0,level], nc[0,level],
+        nprc1tot, npraitot, nprcitot, nsubrtot, nraggtot, nsaggtot, subqr, subnr, \
+        subdr, errstring, prer_evap \
+        = mg.micro_mg_tend(evap_col_steps, evap_steps, col_steps, auto_accr_steps, auto_steps, accr_steps, timestep, t[0,level],
+                           q[0,level], qc[0,level], qi[0,level], nc[0,level],
                            ni[0,level], qr[0,level], qs[0,level], nr[0,level], ns[0,level],
                            relvar_loc[0,level], accre_enhan_loc[0,level], p_loc[0,level], pdel_loc[0,level],
                            precipf_loc[0,level], liqcldf_loc[0,level], icecldf_loc[0,level], naai_loc[0,level],
@@ -386,37 +396,14 @@ def run_mg2(state, level):
     my_ns[0,level] = max(state[ins], 1.e-12)
     tends = mg2_tendencies(level, my_t, my_q, my_qc, my_nc, my_qi, my_ni, my_qr,
                            my_nr, my_qs, my_ns)
-    return s_inv_mult(tends["Total"])
+    output = np.zeros((10*(len(ind)+1),))
+    output[0:10] = s_inv_mult(tends["Total"])
+    for i in range(len(ind)):
+        output[10*(i+1):10*(i+2)] = s_inv_mult(tends[short_names[i]])
+    return output
 
 ind = np.arange(len(short_names))
 
-column = 856
-
-t_loc[:,:] = t[0,:,column]
-q_loc[:,:] = q[0,:,column]
-qc_loc[:,:] = qc[0,:,column]
-qi_loc[:,:] = qi[0,:,column]
-nc_loc[:,:] = nc[0,:,column]
-ni_loc[:,:] = ni[0,:,column]
-qr_loc[:,:] = qr[0,:,column]
-qs_loc[:,:] = qs[0,:,column]
-nr_loc[:,:] = nr[0,:,column]
-ns_loc[:,:] = ns[0,:,column]
-relvar_loc[:,:] = relvar[0,:,column]
-accre_enhan_loc[:,:] = accre_enhan[0,:,column]
-p_loc[:,:] = p[0,:,column]
-pdel_loc[:,:] = pdel[0,:,column]
-precipf_loc[:,:] = precipf[0,:,column]
-liqcldf_loc[:,:] = liqcldf[0,:,column]
-icecldf_loc[:,:] = icecldf[0,:,column]
-naai_loc[:,:] = naai[0,:,column]
-npccn_loc[:,:] = npccn[0,:,column]
-rndst_loc[:,:,:] = rndst[0,:,column:column+1,:].transpose([1, 0, 2])
-nacon_loc[:,:,:] = nacon[0,:,column:column+1,:].transpose([1, 0, 2])
-frzimm_loc[:,:] = frzimm[0,:,column]
-frzcnt_loc[:,:] = frzcnt[0,:,column]
-frzdep_loc[:,:] = frzdep[0,:,column]
-
 def get_modes_at_level(level):
     state = np.zeros((10,))
     state[it] = t_loc[0,level]
@@ -433,15 +420,24 @@ def get_modes_at_level(level):
                            qr_loc, nr_loc, qs_loc, ns_loc)
 
     j_mg2 = ndt.Jacobian(run_mg2, method='forward', order=2)
-    j_loc = s.dot(sq.dot(la.solve_triangular(sr, j_mg2(s_inv_mult(state),
-                                                       level).T, trans='T')).T)
+    j_locs_dimless = j_mg2(s_inv_mult(state), level)
+    j_loc = s.dot(sq.dot(la.solve_triangular(sr, j_locs_dimless[0:10,:].T, trans='T')).T)
     evals, evecs = la.eig(j_loc)
     # Sometimes la.eig spits out complex eigenvalues unnecessarily.
     if all(np.imag(evals) == 0.):
         evals = np.real(evals)
     if (np.imag(evecs) == 0.).all():
         evecs = np.real(evecs)
-    return tends, evals, evecs
+    associations = np.zeros((10, len(ind)))
+    for i in range(len(ind)):
+        j_proc_loc = s.dot(sq.dot(la.solve_triangular(sr, j_locs_dimless[10*(i+1):10*(i+2),:].T, trans='T')).T)
+        j_proc_j_basis = la.solve(evecs, j_proc_loc.dot(evecs))
+        associations[:,i] = np.abs(np.diag(j_proc_j_basis))
+    for i in range(10):
+        assoc_norm = la.norm(associations[i,:], 1)
+        if assoc_norm != 0.:
+            associations[i,:] /= assoc_norm
+    return tends, evals, evecs, associations
 
 def calculate_indices(tends, evals, evecs):
     strengths = la.solve(evecs, tends["Total"])
@@ -490,7 +486,7 @@ def evec_bar_graph(evals, strengths, participations):
 cutoff2 = 0.5
 cutoff3 = 0.9
 tend_cutoff = 1.e-10
-num_columns = 2048
+num_columns = 128
 
 plt.autoscale(tight=True)
 
@@ -514,7 +510,8 @@ def evec_bar_graph(evals, strengths, participations):
 singular_counter = 0
 
 for column in range(num_columns):
-    print("On column: ", column)
+#    if column % 16 == 0:
+    print("On column: ", column, file=sys.stderr)
     t_loc[:,:] = t[0,:,column]
     q_loc[:,:] = q[0,:,column]
     qc_loc[:,:] = qc[0,:,column]
@@ -544,7 +541,7 @@ def evec_bar_graph(evals, strengths, participations):
                                       icecldf_loc, mgncol=1, nlev=lev)
     for level in range(lev):
         c = label[0,level,column]
-        tends, evals, evecs = get_modes_at_level(level)
+        tends, evals, evecs, associations = get_modes_at_level(level)
         tot_t = tends["Total"]
         tend_norm = (abs(tot_t[iq]) + abs(tot_t[iqc]) + abs(tot_t[iqi]) +
                      abs(tot_t[iqr]) + abs(tot_t[iqs])) / 2.
@@ -561,14 +558,14 @@ def evec_bar_graph(evals, strengths, participations):
         for i in range(len(evals)):
             evalues_all[c].append(np.real(evals[i]))
             for j in range(len(short_names)):
-                if np.abs(participations[i,j]) >= cutoff1:
+                if np.abs(associations[i,j]) >= cutoff1:
                     evalues1[c][short_names[j]].append(np.real(evals[i]))
-                    if np.abs(participations[i,j]) >= cutoff2:
+                    if np.abs(associations[i,j]) >= cutoff2:
                         if np.real(evals[i]) < -0.1:
                             evalues2[c][short_names[j]].append(np.real(evals[i]))
                             for j2 in range(len(short_names)):
-                                evalue_correlation[short_names[j]][short_names[j2]] += np.abs(participations[i,j2])
-                        if np.abs(participations[i,j]) >= cutoff3:
+                                evalue_correlation[short_names[j]][short_names[j2]] += np.abs(associations[i,j2])
+                        if np.abs(associations[i,j]) >= cutoff3:
                             evalues3[c][short_names[j]].append(np.real(evals[i]))
 
 #Convert the correlation sum to an average.
@@ -610,7 +607,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.xlabel("Eigenvalue (1/s)")
 plt.xlim(1./max_t, 1./min_t)
 plt.ylabel("Number of modes")
-plt.savefig('./time_hist_all_values_pos.eps')
+plt.savefig('./time_hist_all_values_pos_new.eps')
 plt.close()
 
 # Histogram of all negative eigenvalues.
@@ -631,7 +628,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.xlabel("Eigenvalue (1/s)")
 plt.xlim(1./min_t, 1./max_t)
 plt.ylabel("Number of modes")
-plt.savefig('./time_hist_all_values_neg.eps')
+plt.savefig('./time_hist_all_values_neg_new.eps')
 plt.close()
 
 # Histogram of all near-zero eigenvalues.
@@ -652,7 +649,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.xlim(-1./max_t, 1./max_t)
 plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
 plt.ylabel("Number of modes")
-plt.savefig('./time_hist_all_values_n0.eps')
+plt.savefig('./time_hist_all_values_n0_new.eps')
 plt.close()
 
 # Now broken out by cluster.
@@ -676,7 +673,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.ylabel("Cluster Index")
 plt.clim(vmin=0., vmax=0.2)
 plt.colorbar()
-plt.savefig('./time_hist_cluster_2D_pos.eps')
+plt.savefig('./time_hist_cluster_2D_pos_new.eps')
 plt.close()
 
 # Histogram of all negative eigenvalues.
@@ -697,7 +694,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.ylabel("Cluster Index")
 plt.clim(vmin=0., vmax=0.2)
 plt.colorbar()
-plt.savefig('./time_hist_cluster_2D_neg.eps')
+plt.savefig('./time_hist_cluster_2D_neg_new.eps')
 plt.close()
 
 # Histogram of all near-zero eigenvalues.
@@ -718,7 +715,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.ylabel("Cluster Index")
 plt.clim(vmin=0., vmax=0.2)
 plt.colorbar()
-plt.savefig('./time_hist_cluster_2D_n0.eps')
+plt.savefig('./time_hist_cluster_2D_n0_new.eps')
 plt.close()
 
 # Now broken up by process.
@@ -748,7 +745,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.ylabel("Process")
 plt.clim(vmin=0., vmax=0.2)
 plt.colorbar()
-plt.savefig('./time_hist_process_2D_pos.eps')
+plt.savefig('./time_hist_process_2D_pos_new.eps')
 plt.close()
 
 # Histogram of all negative eigenvalues.
@@ -777,7 +774,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.ylabel("Process")
 plt.clim(vmin=0., vmax=0.2)
 plt.colorbar()
-plt.savefig('./time_hist_process_2D_neg.eps')
+plt.savefig('./time_hist_process_2D_neg_new.eps')
 plt.close()
 
 # Histogram of all near-zero eigenvalues.
@@ -806,7 +803,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.ylabel("Process")
 plt.clim(vmin=0., vmax=0.2)
 plt.colorbar()
-plt.savefig('./time_hist_process_2D_n0.eps')
+plt.savefig('./time_hist_process_2D_n0_new.eps')
 plt.close()
 
 # Correlation array.
@@ -829,7 +826,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.subplots_adjust(left=0.25)
 plt.ylabel("Primary process")
 plt.colorbar()
-plt.savefig('./process_association.eps')
+plt.savefig('./process_association_new.eps')
 plt.close()
 
 
@@ -877,7 +874,7 @@ def evec_bar_graph(evals, strengths, participations):
     ax.tick_params('y', direction='out')
     plt.subplots_adjust(left=0.25)
     plt.colorbar()
-    plt.savefig('./time_hist_2D_pos_c{}.eps'.format(c))
+    plt.savefig('./time_hist_2D_pos_c{}_new.eps'.format(c))
     plt.close()
 
     neg_evalues2 = {}
@@ -912,7 +909,7 @@ def evec_bar_graph(evals, strengths, participations):
     ax.tick_params('y', direction='out')
     plt.subplots_adjust(left=0.25)
     plt.colorbar()
-    plt.savefig('./time_hist_2D_neg_c{}.eps'.format(c))
+    plt.savefig('./time_hist_2D_neg_c{}_new.eps'.format(c))
     plt.close()
 
     n0_evalues2 = {}
@@ -948,7 +945,7 @@ def evec_bar_graph(evals, strengths, participations):
     ax.tick_params('y', direction='out')
     plt.subplots_adjust(left=0.25)
     plt.colorbar()
-    plt.savefig('./time_hist_2D_n0_c{}.eps'.format(c))
+    plt.savefig('./time_hist_2D_n0_c{}_new.eps'.format(c))
     plt.close()
 
 
@@ -963,7 +960,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.title("Number of grid points with given net water mass tendency")
 plt.xlabel("Net water tendency")
 plt.ylabel("Number of grid points")
-plt.savefig('./tend_hist.eps')
+plt.savefig('./tend_hist_new.eps')
 plt.close()
 
 # Log version of plot
@@ -975,7 +972,7 @@ def evec_bar_graph(evals, strengths, participations):
 plt.title("Number of grid points with given net water mass tendency")
 plt.xlabel("Net water tendency")
 plt.ylabel("Number of grid points")
-plt.savefig('./tend_hist_log.eps')
+plt.savefig('./tend_hist_log_new.eps')
 plt.close()
 
 cutoff_score = stats.percentileofscore(tend_norms, tend_cutoff)