Moves NCC diagnostics into separate test file; adds outputs from stru…

…cture to aid diagnostics. Cleans up parts of mdwarf script.

mdwarf_hydro.py

@@ -134,10 +134,6 @@
 lnρ = de.field.Field(dist=d, bases=(b.radial_basis,), dtype=np.float64)
 ρT_inv = de.field.Field(dist=d, bases=(b,), dtype=np.float64)
 
-#T['g'] = structure['T']['g']
-#lnρ['g'] = structure['lnρ']['g']
-logger.info("shape and size of T['g'] {} & {}".format(T['g'].shape, T['g'].size))
-logger.info("shape and size of ρT_inv['g'] {} & {}".format(ρT_inv['g'].shape, ρT_inv['g'].size))
 if T['g'].size > 0 :
     for i, r_i in enumerate(r1[0,0,:]):
          T['g'][:,:,i] = structure['T'](r=r_i).evaluate()['g']
@@ -149,10 +145,7 @@
 ρ = d_exp(lnρ).evaluate()
 grad_lnρ = grad(lnρ).evaluate()
 ρ_inv = d_exp(-lnρ).evaluate()
-ρT_inv_rb = (T_inv*ρ_inv).evaluate()
-ρT_inv_rb.require_scales(1)
-if ρT_inv_rb['g'].size > 0:
-    ρT_inv['g'] = ρT_inv_rb['g']
+ρT_inv = (T_inv*ρ_inv).evaluate()
 
 # Entropy source function, inspired from MESA model
 source = de.field.Field(dist=d, bases=(b,), dtype=np.float64)
@@ -173,49 +166,8 @@
 trace_e = trace(e)
 trace_e.store_last = True
 Phi = trace(dot(e, e)) - 1/3*(trace_e*trace_e)
-# Problem
-if args['--debug']:
-    import matplotlib.pyplot as plt
-    if T['g'].size > 0:
-        fig, ax = plt.subplots(nrows=3, ncols=2)
-        ax[0,0].plot(r[0,0,:], T['g'][0,0,:])
-        ax[0,1].plot(r[0,0,:], lnρ['g'][0,0,:])
-        ax[1,0].plot(r[0,0,:], grad_lnT['g'][2][0,0,:])
-        ax[1,1].plot(r[0,0,:], grad_lnρ['g'][2][0,0,:])
-        ax[2,0].plot(r[0,0,:], T_inv['g'][0,0,:])
-        ax[2,1].plot(r[0,0,:], ρ_inv['g'][0,0,:])
-        ax[0,0].set_ylabel('T')
-        ax[0,1].set_ylabel(r'$\ln \rho$')
-        ax[1,0].set_ylabel('gradT')
-        ax[1,1].set_ylabel('gradlnrho')
-        ax[2,0].set_ylabel('1/T')
-        ax[2,1].set_ylabel('1/ρ')
-        plt.tight_layout()
-        fig.savefig('nccs_p{}.pdf'.format(rank))
-        print(grad_lnρ['g'][2])
-        print("rho inv: {}".format(ρ_inv['g']))
-
-    if T['c'].size > 0:
-        fig, ax = plt.subplots(nrows=3, ncols=2)
-        ax[0,0].plot(np.abs(T['c'][0,0,:]))
-        ax[0,1].plot(np.abs(lnρ['c'][0,0,:]))
-        ax[1,0].plot(np.abs(grad_lnT['c'][1][0,0,:])) # index 1 is spin 0
-        ax[1,1].plot(np.abs(grad_lnρ['c'][1][0,0,:])) # index 1 is spin 0
-        ax[2,0].plot(np.abs(T_inv['c'][0,0,:]))
-        ax[2,1].plot(np.abs(ρ_inv['c'][0,0,:]))
-        ax[0,0].set_ylabel('T')
-        ax[0,1].set_ylabel(r'$\ln \rho$')
-        ax[1,0].set_ylabel('gradT')
-        ax[1,1].set_ylabel('gradlnrho')
-        ax[2,0].set_ylabel('1/T')
-        ax[2,1].set_ylabel('1/ρ')
-        for axi in ax:
-            for axii in axi:
-                axii.axhline(y=ncc_cutoff, linestyle='dashed', color='xkcd:grey')
-                axii.set_yscale('log')
-        plt.tight_layout()
-        fig.savefig('nccs_coeff_p{}.pdf'.format(rank))
 
+#Problem
 problem = problems.IVP([u, p, s, τ_u, τ_s], ncc_cutoff=ncc_cutoff)
 problem.add_equation((ddt(u) + grad(p) - Co2*T*grad(s) - Ek*ρ_inv*viscous_terms + LiftTau(τ_u,-1),
                       - dot(u, e) - cross(ez_g, u)), condition = "ntheta != 0")
@@ -273,10 +225,10 @@
         T0 = np.sum(vol_correction*weight_r*weight_theta*0.5*s['g']**2)
         T0 *= (np.pi)/(Lmax+1)/L_dealias
         T0 = reducer.reduce_scalar(T0, MPI.SUM)
-        logger.info("iter: {:d}, dt={:e}, t={:e}, E0={:e}, T0={:e}".format(solver.iteration, dt, solver.sim_time, E0, T0))
+        logger.info("iter: {:d}, dt={:.2e}, t={:.3e}, E0={:e}, T0={:e}".format(solver.iteration, dt, solver.sim_time, E0, T0))
         good_solution = np.isfinite(E0)
     elif solver.iteration % report_cadence == 0:
-        logger.info("iter: {:d}, dt={:e}, t={:e}".format(solver.iteration, dt, solver.sim_time))
+        logger.info("iter: {:d}, dt={:.2e}, t={:.3e}".format(solver.iteration, dt, solver.sim_time))
     if solver.iteration % hermitian_cadence in timestepper_history:
         for field in solver.state:
             field['g']

structure.py

@@ -47,6 +47,8 @@ def error(perts):
     structure = {'T':T,'lnρ':lnρ}
     for key in structure:
         structure[key].require_scales(1)
+    structure['r'] = r
+    structure['problem'] = {'c':c, 'b':b, 'problem':problem}
     return structure
 
 if __name__=="__main__":

test_ncc_scaling.py

@@ -0,0 +1,96 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from structure import lane_emden
+import dedalus.public as de
+
+import logging
+logger = logging.getLogger(__name__)
+dlog = logging.getLogger('matplotlib')
+dlog.setLevel(logging.WARNING)
+
+ncc_cutoff = 1e-10
+
+LE = lane_emden(63)
+
+T = LE['T']
+lnρ = LE['lnρ']
+c = LE['problem']['c']
+r = LE['r']
+
+d_exp = lambda A: de.operators.UnaryGridFunction(np.exp, A)
+d_log = lambda A: de.operators.UnaryGridFunction(np.log, A)
+power = lambda A, B: de.operators.Power(A, B)
+grad = lambda A: de.operators.Gradient(A, c)
+
+lnT = d_log(T).evaluate()
+T_inv = power(T,-1).evaluate()
+grad_lnT = grad(lnT).evaluate()
+ρ = d_exp(lnρ).evaluate()
+grad_lnρ = grad(lnρ).evaluate()
+ρ_inv = d_exp(-lnρ).evaluate()
+ρT_inv = (T_inv*ρ_inv).evaluate()
+
+# the LHS anelastic NCCs are T, grad_lnT, grad_lnρ, ρ_inv
+
+fig, ax = plt.subplots(nrows=2, ncols=2)
+ax[0,0].plot(r[0,0,:], T['g'][0,0,:])
+ax[0,1].plot(r[0,0,:], ρ_inv['g'][0,0,:])
+ax[1,0].plot(r[0,0,:], grad_lnT['g'][2][0,0,:])
+ax[1,1].plot(r[0,0,:], grad_lnρ['g'][2][0,0,:])
+ax[0,0].set_ylabel('T')
+ax[0,1].set_ylabel('1/ρ')
+ax[1,0].set_ylabel('gradT')
+ax[1,1].set_ylabel('gradlnrho')
+plt.tight_layout()
+fig.savefig('nccs.pdf')
+
+fig, ax = plt.subplots(nrows=2, ncols=2)
+ax[0,0].plot(np.abs(T['c'][0,0,:]))
+ax[0,1].plot(np.abs(ρ_inv['c'][0,0,:]))
+ax[1,0].plot(np.abs(grad_lnT['c'][1][0,0,:])) # index 1 is spin 0
+ax[1,1].plot(np.abs(grad_lnρ['c'][1][0,0,:])) # index 1 is spin 0
+ax[0,0].set_ylabel('T')
+ax[0,1].set_ylabel('1/ρ')
+ax[1,0].set_ylabel('gradT')
+ax[1,1].set_ylabel('gradlnrho')
+for axi in ax:
+    for axii in axi:
+        axii.axhline(y=ncc_cutoff, linestyle='dashed', color='xkcd:grey')
+        axii.set_yscale('log')
+plt.tight_layout()
+fig.savefig('nccs_coeff.pdf')
+
+
+
+T_1 = (1/T).evaluate()
+grad_lnT_1 = (1/grad_lnT).evaluate()
+grad_lnρ_1 = (1/grad_lnρ).evaluate()
+ρ_inv_1 = (1/ρ_inv).evaluate()
+
+fig, ax = plt.subplots(nrows=2, ncols=2)
+ax[0,0].plot(r[0,0,:], T_1['g'][0,0,:])
+ax[0,1].plot(r[0,0,:], ρ_inv_1['g'][0,0,:])
+ax[1,0].plot(r[0,0,:], grad_lnT_1['g'][2][0,0,:])
+ax[1,1].plot(r[0,0,:], grad_lnρ_1['g'][2][0,0,:])
+ax[0,0].set_ylabel('T')
+ax[0,1].set_ylabel('1/ρ')
+ax[1,0].set_ylabel('gradT')
+ax[1,1].set_ylabel('gradlnrho')
+plt.tight_layout()
+fig.savefig('nccs_1.pdf')
+
+fig, ax = plt.subplots(nrows=2, ncols=2)
+ax[0,0].plot(np.abs(T_1['c'][0,0,:]))
+ax[0,1].plot(np.abs(ρ_inv_1['c'][0,0,:]))
+ax[1,0].plot(np.abs(grad_lnT_1['c'][1][0,0,:])) # index 1 is spin 0
+ax[1,1].plot(np.abs(grad_lnρ_1['c'][1][0,0,:])) # index 1 is spin 0
+ax[0,0].set_ylabel('T')
+ax[0,1].set_ylabel('1/ρ')
+ax[1,0].set_ylabel('gradT')
+ax[1,1].set_ylabel('gradlnrho')
+for axi in ax:
+    for axii in axi:
+        axii.axhline(y=ncc_cutoff, linestyle='dashed', color='xkcd:grey')
+        axii.set_yscale('log')
+plt.tight_layout()
+fig.savefig('nccs_1_coeff.pdf')