Changes between Version 81 and Version 82 of FluxLimitedDiffusion

Timestamp:

03/27/13 22:07:16 (12 years ago)

Author:

Jonathan

Comment:

—

FluxLimitedDiffusion

@@ -89 +89 @@
 For static diffusion, the terms in blue with v^2^/c can be dropped and the system can be split into the usual hydro update (black), radiative source terms (green), and a coupled implicit solve (red) for the radiation energy density and thermal energy density (ie temperature). 
 
-[[CollapsibleStart(Operator Splitting)]]
+[[CollapsibleStart(Operator Splitting 1)]]
 == Operator splitting ==
 Krumholz et al. perform Implicit Radiative, Explicit Hydro, Explicit Radiative
@@ -135 +135 @@
 [[CollapsibleEnd()]]
 
-[[CollapsibleStart(Explcit Update)]]
+[[CollapsibleStart(Explcit Update 1)]]
+== Eplicit Update 1 ==
 The extra terms in the explicit update due to radiation energy are as follows:
 
@@ -174 +175 @@
 [[CollapsibleEnd()]]
 
-[[CollapsibleStart(Implicit Update)]]
-== Implicit Equation ==
+[[CollapsibleStart(Implicit Update 1)]]
+== Implicit Update 1 ==
 
 For now we will assume that [[latex(\kappa_{0P})]] and [[latex(\kappa_{0R})]] are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions.  In this case we can solve the radiation energy equations:
@@ -210 +211 @@
 
 
-=== Implicit Discretization ===
+=== Implicit Discretization 1 ===
 Which we can discretize for (1D) as
 
@@ -339 +340 @@
 
 [[CollapsibleEnd()]]
+
+  [[latex(\frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P\color{green}{-\lambda \nabla E})]]   
+  [[latex(\frac{\partial e}{\partial t}  + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = \color{red}{-\kappa_{0P}(4 \pi B-cE)} \color{orange}{+\lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E} \color{blue}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E})]]  
+  [[latex(\frac{\partial E}{\partial t}  \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E} = \color{red}{\kappa_{0P} (4 \pi B-cE)} \color{orange}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )} \color{blue}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}  )]]  
+
+For static diffusion, the terms in blue with v^2^/c can be dropped and the system can be split into the usual hydro update (black), radiative source terms (green), and a coupled implicit solve (red) for the radiation energy density and thermal energy density (ie temperature). And the terms in orange can be solved semi-implicitly (velocity terms are explicit, while radiation energy terms are implicit) 
+
+
+[[CollapsibleStart(Operator Splitting 2)]]
+== Operator splitting 2 ==
+Krumholz et al. perform Implicit Radiative, Explicit Hydro, Explicit Radiative
+
+In AstroBEAR this would look like:
+* Initialization
+ * Prolongate, d, p, e, E, Edot
+* Step 1
+ * Overlap d, p, e, E and do physical BC's
+ * Do IR which updates e,,0,,, and E,,0,, using d,,1,,, e,,1,,, E,,1,,
+ * Update E,,2*mbc,, using Edot,,2*mbc,,
+ * Update e,,2*mbc,, using E,,2*mbc,,, Edot,,2*mbc,,, and e,,2*mbc,,
+ * Update Edot,,0,, using pre IR and post IR E,,0,,
+ * Ghost e,,2*mbc,,, E,,mbc+1,,, Edot,,mbc+1,,
+ * Do first EH,,mbc,,
+ * Do ER,,mbc,, --- Terms with grad E can be done without ghosting since EH did not change E.  The del dot vE term needs time centered face centered velocities which can be stored during the hydro update.
+ * Store Edot in child arrays to be prolongated
+* Step 2
+ * Overlap d, p, e, E, and do physical BCs
+ * Do IR which updates e,,0,,, and E,,0,, using d,,1,,, e,,1,,, E,,1,,
+ * Update Edot,,0,, using pre IR and post IR E,,0,,
+ * Update E,,1,, using Edot,,1,,
+ * Ghost e,,mbc,,, E,,1,,, Edot,,1,,
+ * Do second EH,,0,,
+ * Do ER,,0,, --- Terms with grad E can be done without ghosting since EH did not change E.  The del dot vE term needs time centered face centered velocities which can be stored during the hydro update.
+
+
+However with AMR is we need coarse boundary values for E at t=0, and t=dt, and t = 2dt and we would like coarse boundary values at t=2dt to match solution from coarse grid...
+
+But the coarse update involves one implicit and one explicit solve.  Is there a way to interpolate the fine grid ghost zones in time without storing two time derivatives?  If we store dE_i+dE_e then we could do an explicit radiative update internally (E^n^ -> E*), then an implicit radiative update (E* -> E^n+1^ using the fully updated ghost zones (E^n+1^,,g,, = E^n^,,g,,+dE,,g,,) and then store the new time derivative (dE=E^n+1^-E^n^)
+
+This seems like it would work for the radiative field, but what about the internal energy terms?  We have an equation for e^n+1^ that is a function of E^n+1^, e* and E* but we have no way of getting E* in the ghost zones...
+
+So the solution is perhaps to store the contributions from the implicit update and the non-conservative heating terms that appear in the internal energy equation.  Then we can update E in the ghost zones for the implicit and the first explicit term using the time derivative - and then use the new E to update e... Then the conservative explicit RadEnergy term can be calculated after the hydro step as well as the momentum explicit term - and then this flux can be coarsened - to keep the value of E in the coarse cells consistent with the value of E in the fine cells...  Since both will have been updated by the same eDot and fluxes...
+
+
+[[CollapsibleEnd()]]
+
+[[CollapsibleStart(Explcit Update 2)]]
+== Eplicit Update 2 ==
+The extra terms in the explicit update due to radiation energy are as follows:
+
+ [[latex(\frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) =-\lambda \nabla E)]]   
+ [[latex(\frac{\partial e}{\partial t}  = \lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E )]]  
+ [[latex(\frac{\partial E}{\partial t}  = -\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E -\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right ) )]]  
+
+These can be discretized as follows:
+
+ [[latex(p^{n+1}_i=p^n_i - \frac{1}{2}\frac{\Delta t}{\Delta x} \lambda_{i} \left ( E^n_{i+1}-E^n_{i-1} \right ) )]]
+
+ [[latex(e^{n+1}_i=e^n_i + \frac{1}{2}\frac{\Delta t}{\Delta x} \lambda_i \left ( 2 \frac{\kappa^n_{i,0P}}{\kappa^n_{i,0R}}-1 \right ) \left ( v^n_i \left ( E^n_{i+1}-E^n_{i-1} \right ) \right ) )]]  
+
+ [[latex(E^{n+1}_i=E^n_i - \frac{\Delta t}{\Delta x} \left ( \frac{\lambda_i}{2} \left ( 2 \frac{\kappa^n_{i,0P}}{\kappa^n_{i,0R}}-1 \right ) \left ( v^n_i \left ( E^n_{i+1}-E^n_{i-1} \right ) \right ) + \left ( F_{i+1/2}-F_{i-1/2} \right ) \right ))]]  
+
+where
+[[latex(F_{i+1/2}=\frac{3-R_{2,i+1/2}}{8} \left(v_{i}^n+v_{i+1}^n \right ) \left ( E^n_i+E^n_{i+1} \right ) )]]
+
+where
+
+[[latex(R_{2,i+1/2} = \lambda_{i+1/2}+\lambda_{i+1/2}^2 R_{i+1/2}^2)]]
+and 
+[[latex(R_{i+1/2} = \frac{\left | E^n_{i+1}-E^n_{i} \right | }{2 \kappa_{0R,i+1/2} \left ( E^n_i+E^n_{i+1} \right )})]]
+
+and
+[[latex(\lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right ) )]]
+
+and 
+[[latex(\kappa_{0R,i+1/2} = \frac{\kappa^n_{0R,i}+\kappa^n_{0R,i+1}}{2} \mbox{ and } \lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right ) )]]
+
+and 
+
+[[latex(\lambda_{i} = \frac{1}{R_{i}} \left ( \coth R_{i} - \frac{1}{R_{i}} \right ) )]]
+
+and
+[[latex(R_{i} = \frac{\left | E^n_{i+1}-E^n_{i-1} \right | }{2 \kappa_{0R,i} E^n_{i}})]]
+
+[[CollapsibleEnd()]]
+
+[[CollapsibleStart(Implicit Update 1)]]
+== Implicit Update 1 ==
+
+For now we will assume that [[latex(\kappa_{0P})]] and [[latex(\kappa_{0R})]] are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions.  In this case we can solve the radiation energy equations:
+
+   [[latex(\frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B(T)-cE) = \nabla \cdot \mathbf{F} + \kappa_{0P} (4 \pi B(T)-cE))]]   
+   [[latex(\frac{\partial e}{\partial t} = - \kappa_{0P} (4 \pi B(T)-cE))]]   
+
+where [[latex(\mathbf{F} = \frac{c\lambda}{\kappa_{0R}} \nabla E)]]
+
+=== Expanding about e,,0,, ===
+
+Of course even if the opacity is independent of energy and radiation energy, the above combined system of equations is still non-linear due to the dependence on Temperature of the Planck Function [[latex(B(T))]]
+
+If we ignore the changes in the Temperature due to heating during the implicit step which would feed back into the source function.  We can improve this by writing 
+
+[[latex(B(T) = B \left ( T_0+dT \right ) = B \left ( T_0 \right ) + \left . \frac{\partial B}{\partial T} \right | _{T_0} \frac{\partial T}{\partial e} de = \frac{c}{4 \pi} a_R \left ( T_0^4 + 4T_0^3\Gamma de \right ) = B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right ) )]]
+
+where
+
+[[latex(\Gamma = \frac{\partial T}{\partial e} = \frac{(\gamma-1)}{n k_B})]]
+
+Then the system of equations becomes
+
+   [[latex(\frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} \left [4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] )]]   
+   [[latex(\frac{\partial e}{\partial t} = - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ] )]]   
+
+which will be accurate as long as [[latex(4\Gamma \frac{e-e_0}{T_0} < \xi << 1)]] or [[latex(\Delta e = e-e_0 < \xi \frac{T_0}{4 \Gamma})]]
+
+We can calculate the time scale for this to be true using the evolution equation for the energy density
+
+[[latex(\Delta e = -\Delta t \kappa_{0P} \left [ 4 \pi B_0 -cE \right ] < \xi \frac{T_0}{4 \Gamma})]]
+
+which gives [[latex(\Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left ( 4 \pi B_0 - cE \right ) })]]
+
+
+=== Implicit Discretization 2 ===
+Which we can discretize for (1D) as
+
+   [[latex(E^{n+1}_i-E^{n}_i = \left [ \alpha^n_{i+1/2} \left ( E^{*}_{i+1}-E^{*}_{i} \right ) - \alpha^n_{i-1/2} \left ( E^{*}_{i}-E^{*}_{i-1} \right ) \right ] - \epsilon^n_i E^{*}_i  + \phi^n_i e^{*}_i  + \theta^n_i) )]]   
+   [[latex(e^{n+1}_i-e^{n}_i = \epsilon^n_i E^{*}_i  - \phi^n_i e^{*}_i  - \theta^n_i )]]   
+
+where the diffusion coefficient is given by 
+
+[[latex(\alpha_{i+1/2}=\frac{\Delta t}{\Delta x^2}  \frac{c \lambda_{i+1/2}}{\kappa_{0R,i+1/2}} \mbox{ where } \kappa_{0R,i+1/2} = \frac{\kappa^n_{0R,i}+\kappa^n_{0R,i+1}}{2} \mbox{ and } \lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right ) )]]
+and where
+[[latex(R_{i+1/2} = \frac{\left | E^n_{i+1}-E^n_{i} \right | }{2 \kappa_{0R,i+1/2} \left ( E^n_i+E^n_{i+1} \right )})]]
+
+
+and 
+
+[[latex(\epsilon^n_i=c\Delta t \kappa^n_{0P,i})]]
+
+represents the number of absorption/emissions during the time step
+
+
+and 
+
+[[latex(\phi = \epsilon^n_i \frac{4 \pi}{c} B \left ( T^n_i \right ) \left ( \frac{4\Gamma}{T^n_i} \right ) )]]
+
+[[latex(\theta = \epsilon^n_i \frac{4 \pi}{c} B \left ( T^n_i \right ) \left ( 1 - 4\Gamma \frac{e^n_i}{T^n_i} \right ) )]]
+
+
+and we can think of the radiative flux as
+
+[[latex(\frac{\Delta t}{\Delta x}\mathbf{F}^n_{i+1/2} = \alpha^n_{i+1/2} \left ( E^{*}_{i+1} - E^{*}_i \right ) )]]
+
+=== Time Discretization ===
+
+Now all the terms on the right hand side that are linear in E or e have been written as E^*^ or e^*^ because there are different ways to approximate E^*^ (e^*^).  For Backward Euler we have
+[[latex(E^*_i = E^{n+1}_i)]]
+and for Crank Nicholson we have
+[[latex(E^*_i = \frac{1}{2} \left ( E^{n+1}_i + E^n_i \right ) )]]
+or we can parameterize the solution
+[[latex(E^*_i = \psi E^{n+1}_i + \bar{\psi}E^n_i)]]
+where [[latex(\bar{\psi} = 1-\psi)]]
+
+Backward Euler has [[latex(\psi=1)]] and Crank Nicholson has [[latex(\psi=1/2)]]
+
+Forward Euler has [[latex(\psi=0)]]
+
+In any event in 1D we have the following matrix coefficients
+
+   [[latex(\left [ 1 + \psi \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^{n+1}_i - \left ( \psi \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1} - \left ( \psi \phi^n_i \right ) e^{n+1}_i=\left [ 1 - \bar{\psi} \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha^n_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha^n_{i-1/2} \right ) E^{n}_{i-1} +\bar{\psi}\phi^n_i e^n_i + \theta^n_i)]]   
+   [[latex(\left ( 1 +\psi \phi^n_i \right ) e^{n+1}_i - \left ( \psi \epsilon^n_i \right )E^{n+1}_i =\left ( 1 - \bar{\psi}\phi^n_i \right ) e^n_i + \left ( \bar{\psi} \epsilon^n_i \right ) E^n_i-\theta^i_n )]]   
+
+
+Now since the second equation has no spatial dependence, we can solve it for 
+   [[latex(\color{purple}{e^{n+1}_i = \frac{1}{ 1 +\psi \phi^n_i}\left \{ \left ( \psi \epsilon^n_i \right )E^{n+1}_i + \left ( 1 - \bar{\psi}\phi^n_i \right ) e^n_i + \left ( \bar{\psi} \epsilon^n_i \right ) E^n_i-\theta^i_n \right \}} )]]   
+
+and plug the result into the first equation to get a matrix equation involving only one variable.
+
+   [[latex(\color{purple}{\left [ 1 + \psi \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \frac{\epsilon^n_i}{ 1 +\psi \phi^n_i}\right ) \right ] E^{n+1}_i - \left ( \psi \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2}  +\frac{\epsilon^n_i }{ 1 +\psi \phi^n_i} \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha^n_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha^n_{i-1/2} \right ) E^{n}_{i-1} + \frac{ \phi^n_i}{ 1 +\psi \phi^n_i}  e^n_i+ \frac{1}{ 1 +\psi \phi^n_i}\theta^i_n})]]   
+
+
+
+If we ignore the change in the Planck function due to heating during the implicit solve, it is equivalent to setting [[latex(\psi \phi = 0)]]  This gives the following equations:
+
+   [[latex(\left [ 1 + \psi \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^{n+1}_i - \left ( \psi \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^n_i+ \left ( \bar{\psi} \alpha^n_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha^n_{i-1/2} \right ) E^{n}_{i-1} +\phi^n_i e^n_i + \theta^n_i)]]   
+   [[latex(e^{n+1}_i = e^n_i + \epsilon^n_i  \left [ \left ( \psi E^{n+1}_i + \bar{\psi} E^{n}_i \right ) - \frac{4 \pi}{c} B \left ( T^n_i \right )  \right ] )]]   
+
+In this case the first equation decouples and can be solved on it's own, and then the solution plugged back into the second equation to solve for the new energy.  
+
+
+=== 2D etc... ===
+
+For 2D or 3D we have more connections to add to the matrix elements but it is very straight forward...  There will be additional alpha terms for each  dimension, but everything else stays the same.
+
+
+=== Initial solution vector ===
+For the initial solution vector, we can just use Edot from the parent update (or last time step if we are on the coarse grid) to guess E, and then we can solve for the new e given our guess at the new E using the same time stepping (Backward Euler, Crank Nicholson, etc...).
+
+[[latex(E^{n+1}_i = E^{n}_i+\dot{E}^{n}_i \Delta t)]]
+[[latex(e^{n+1}_i = \frac{1}{1+\psi \phi^n_i} \left \{ \left ( \psi \epsilon^n_i \right )E^{n+1}_i + \left ( 1 - \bar{\psi}\phi^n_i \right ) e^n_i + \left ( \bar{\psi} \epsilon^n_i \right ) E^n_i-\theta^i_n \right \} )]]
+
+== Coarse-Fine Boundaries ==
+ Since we are doing our implicit solves first, we can use time interpolated solutions for the implicit solve for non-refined ghost zones.  To do this we just need Edot.  The opacities etc... in the ghost zones can be obtained from the hydro terms.  
+
+ And the radiative implicit heating in coarse ghost cells can be done along with the initial solution vector so they are available for the hydro update.
+
+
+== Physical Boundary Conditions ==
+
+=== Open (Free streaming) boundaries ===
+We would like the radiation to leave at the free streaming limit.  So 
+[[latex(\frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x})]]
+
+Clearly if we set [[latex(E_g = 0)]] and [[latex(\lambda_g =\kappa_g \Delta x)]] we should get the correct flux.
+
+This corresponds to an
+[[latex(\alpha = c \frac{\Delta t}{\Delta x})]]
+
+So we would just modify [[latex(\alpha)]] and zero out the matrix coefficient to the ghost zone
+
+=== Constant Slope Boundary ===
+Here we want the flux to be constant so energy does not pile up near the boundary.  If we cancel all derivative terms on both sides of the cell, this will effectively match the incoming flux with the outgoing flux.  This can also be done by setting [[latex(\alpha_g = \alpha_i = 0)]]
+
+=== Periodic Boundary ===
+This is the same as internal zones - it just maps the neighbor cell to be across the domain.  Hypre has built in functionality for this under for the Struct Interface
+
+=== User defined radiation field/Coarse Fine boundary ===
+
+This will be the boundary at internal coarse-fine boundaries, but could also be used at the physical boundary if the radiation energy were specified.
+
+=== Reflecting/ZeroSlope Boundary ===
+Reflecting boundary should be fairly straightforward.  This an be achieved by setting [[latex(\alpha_g = 0)]] which zeros out any flux - and has the same effect as setting [[latex(E^{*}_g=E^{*}_i)]] or [[latex(E^{n+1}_g=E^{n+1}_i \mbox{ and } E^{n}_g=E^{n}_i)]]
+
+=== Constant radiative flux ===
+To have a constant radiative flux we must zero out terms involving the gradient and just add [[latex(F_0 \frac{\Delta t}{\Delta x})]] in the source vector
+
+=== Summary ===
+
+
+|| Numerical value  ||  Boundary                  ||  [[latex(E^{n+1}_{i+1})]]        ||  [[latex(E^{n+1}_i)]]                         ||  [[latex(S)]] 
+||  0               ||  RAD_FREE_STREAMING        ||  [[latex(0)]]                    ||  [[latex(\psi c \frac{\Delta t}{\Delta x})]]  ||  [[latex(-\bar{\psi}c \frac{\Delta t}{\Delta x} E^n_i)]]    || 
+||  1               ||  RAD_EXTRAPOLATE_FLUX*     ||  [[latex(0)]]                    ||  [[latex(0)]]                                 ||  [[latex(0)]]    ||
+||  2               ||  INTERNAL/PERIODIC         ||  [[latex(-\psi\alpha_{i+1/2})]]  ||  [[latex(\psi\alpha_{i+1/2})]]                ||  [[latex(\bar{\psi}\alpha_{i+1/2} \left ( E^n_{i+1}-E^n_i \right ))]]  ||
+||  3               ||  RAD_REFLECTING            ||  [[latex(0)]]                    ||  [[latex(0)]]                                 ||  [[latex(0)]]    ||
+||  4               ||  RAD_USERDEFINED_BOUNDARY/AMR_BOUNDARY  ||  [[latex(0)]]       ||  [[latex(\psi\alpha_{i+1/2})]]                ||   [[latex(\bar{\psi}\alpha_{i+1/2} \left ( E^n_{i+1}-E^n_i \right ) - \psi \alpha_{i+1/2} E^{n+1}_{i+1} )]]   ||
+||  5               ||  RAD_USERDEFINED_FLUX      ||  [[latex(0)]]                    ||  [[latex(0)]]                                 ||     [[latex(F_0 \frac{\Delta t}{\Delta x})]]   ||
+
+ !* Also requires zeroing out stencil and source contributions from [[latex(E^{n+1}_{i-1} \mbox{ and } E^n_{i-1})]] 
+
+
+[[CollapsibleEnd()]]