Montecarlo/fix convergence

docs/montecarlo.rst

@@ -10,4 +10,56 @@ Different line interactions
 
 line_interaction_id == 0: scatter
 line_interaction_id == 1: downbranch
-line_interaction_id == 2: macro
+line_interaction_id == 2: macro
+
+Radiationfield estimators
+-------------------------
+
+During the monte-carlo run we collect two estimators for the radiation field:
+
+.. math::
+
+    J_\textrm{estimator} &= \sum{\epsilon l}\\
+    \bar{\nu}_\textrm{estimator} &=  \sum{\epsilon \nu l},
+
+where :math:`\epsilon, \nu` are comoving energy and comoving frequency of a packet respectively.
+
+To calculate the temperature and dilution factor we first calculate the mean intensity in each cell
+( :math:`J = \frac{1}{4\pi\, \Delta t\, V} J_\textrm{estimator}` )., :cite:`2003A&A...403..261L`.
+
+The weighted mean frequency is used to obtain the radiation temperature. Specifically, the radiation temperature is chosen as the 
+temperature of a black body that has the same weighted mean frequency as has been computed in the simulation. Accordingly,
+
+.. math::
+
+    \frac{h \bar{\nu}}{k_{B} T_{R}} = \frac{h}{k_{B} T_{R}} \frac{\bar{\nu}_\textrm{estimator}}{J_\textrm{estimator}} 
+      = 24 \zeta(5) \frac{15}{\pi^4},
+
+where the evaluation comes from the mean value of
+
+.. math::
+
+    \bar{x} = \frac{ \int_0^{\infty} x^4 dx/ (\exp{x} - 1)}{\int_0^{\infty} x^3 dx/ (\exp{x} - 1)} = \frac{24 \zeta(5)}{\pi^4/15}
+
+and so
+
+.. math::
+
+    T_{R} &= \frac{\pi^4}{15} \frac{1}{24 \zeta(5)} \frac{h}{k_{B}} \frac{\bar{\nu}_\textrm{estimator}}{J_\textrm{estimator}} \\
+    &= 0.260945 \frac{h}{k_{B}} \frac{\bar{\nu}_\textrm{estimator}}{J_\textrm{estimator}}.
+
+With the radiation temperature known, we can then obtain our estimate for for the dilution factor. Our radiation field model in the 
+nebular approximation is
+
+.. math::
+
+    J = W B(T_{R}) = W \frac{\sigma_{SB}}{\pi} T_{R}^4,
+
+i.e. a dilute blackbody. Therefore we use our value of the mean intensity derrived from the estimator (above) to obtain the 
+dilution factor
+
+.. math::
+
+    W = \frac{\pi J}{\sigma_{SB} T_{R}^4} = \frac{1}{4\sigma_{SB} T_{R}^4\, \Delta t\, V} J_\textrm{estimator}.
+
+There endeth the lesson.

tardis/config_reader.py

@@ -179,7 +179,6 @@ def luminosity_outer(self):
     @luminosity_outer.setter
     def luminosity_outer(self, value):
         self._luminosity_outer = value
-        self.time_of_simulation = 1 / self.luminosity_outer
 
     def set_densities(self, densities):
         """

tardis/model_radial_oned.py

@@ -13,7 +13,7 @@
 h = constants.cgs.h.value
 kb = constants.cgs.k_B.value
 
-trad_estimator_constant = 0.260944706 * h / kb # (pi**4 / 15)/ (24*zeta(5))
+
 
 w_estimator_constant = (c ** 2 / (2 * h)) * (15 / np.pi ** 4) * (h / kb) ** 4 / (4 * np.pi)
 
@@ -118,7 +118,6 @@ def __init__(self, configuration_object, atom_data):
         self.create_packets()
 
 
-
     #Selecting plasma class
         self.plasma_type = configuration_object.plasma_type
         if self.plasma_type == 'lte':
@@ -192,6 +191,15 @@ def line_interaction_type(self, value):
         self.atom_data.prepare_atom_data(self.selected_atomic_numbers,
             line_interaction_type=self.line_interaction_type, max_ion_number=None)
 
+    @property
+    def t_inner(self):
+        return self._t_inner
+
+    @t_inner.setter
+    def t_inner(self, value):
+        self._t_inner = value
+        self.luminosity_inner = 4 * np.pi * constants.cgs.sigma_sb.value * self.r_inner[0]**2 * self._t_inner**4
+        self.time_of_simulation = 1 / self.luminosity_inner
 
 
     def create_packets(self):
@@ -241,11 +249,34 @@ def initialize_plasmas(self):
 
             # update plasmas
 
-    def calculate_updated_trads(self, nubar_estimators, j_estimators):
-        return trad_estimator_constant * nubar_estimators / j_estimators
+    def calculate_updated_radiationfield(self, nubar_estimator, j_estimator):
+        """
+        Calculate an updated radiation field from the :math:`\\bar{nu}_\\textrm{estimator}` and :math:`\\J_\\textrm{estimator}`
+        calculated in the montecarlo simulation. The details of the calculation can be found in the documentation.
+
+        Parameters
+        ----------
+
+        nubar_estimator : ~np.ndarray (float)
+
+        j_estimator : ~np.ndarray (float)
+
+        Returns
+        -------
+
+        updated_t_rads : ~np.ndarray (float)
+
+        updated_ws : ~np.ndarray (float)
+
+        """
+        trad_estimator_constant = 0.260944706 * h / kb # (pi**4 / 15)/ (24*zeta(5))
+
+        updated_t_rads = trad_estimator_constant * nubar_estimator / j_estimator
+        updated_ws = j_estimator / (4 * constants.cgs.sigma_sb.value * updated_t_rads**4 * self.time_of_simulation * self.volumes)
+
+        return updated_t_rads, updated_ws
+
 
-    def calculate_updated_ws(self, j_estimators, updated_t_rads):
-        return w_estimator_constant * j_estimators / (self.time_of_simulation * (updated_t_rads ** 4) * self.volumes)
 
 
     def update_plasmas(self, updated_t_rads, updated_ws=None):

tardis/simulation.py

@@ -34,8 +34,7 @@ def run_single_radial1d(radial1d_model):
 def run_radial1d(radial1d_model):
     for i in range(9):
         out_nu, out_energy, j_estimators, nubar_estimators =  montecarlo_multizone.montecarlo_radial1d(radial1d_model)
-        updated_t_rads = radial1d_model.calculate_updated_trads(nubar_estimators, j_estimators)
-        updated_ws = radial1d_model.calculate_updated_ws(j_estimators, updated_t_rads)
+        updated_t_rads, updated_ws = radial1d_model.calculate_updated_radiationfield(nubar_estimators, j_estimators)
 
 
         new_trads = 0.5 * (radial1d_model.t_rads + updated_t_rads)
@@ -46,8 +45,10 @@ def run_radial1d(radial1d_model):
 
 
         emitted_energy = radial1d_model.emitted_inner_energy * np.sum(out_energy[out_energy >= 0]) / 1.
+        absorbed_energy = radial1d_model.emitted_inner_energy * np.sum(out_energy[out_energy < 0]) / 1.
 
-        print "Luminosity emitted = %.5e Luminosity requested = %.5e" % (emitted_energy, radial1d_model.luminosity_outer)
+        print "Luminosity emitted = %.5e Luminosity absorbed = %.5e Luminosity requested = %.5e" % (emitted_energy,
+                                                                                    absorbed_energy,radial1d_model.luminosity_outer)
         new_t_inner = radial1d_model.t_inner * (emitted_energy / radial1d_model.luminosity_outer)**-.25
         print "new t_inner = %.2f" % (new_t_inner,)