Python Mesh.getCellCenters方法代码示例

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
          
          # run the rad-hydro transient
          rad_new, hydro_new = runNonlinearTransient(
             mesh         = mesh,
             problem_type = 'rad_hydro',
             dt_option    = 'constant',
             dt_constant  = dt_value,
             slope_limiter = limiter,
             time_stepper = 'BDF2',
             use_2_cycles = True,
             t_start      = t_start,
             t_end        = t_end,
             rad_BC       = rad_BC,
             cross_sects  = cross_sects,
             rad_IC       = rad_IC,
             hydro_IC     = hydro_IC,
             hydro_BC     = hydro_BC,
             mom_src      = mom_src,
             E_src        = E_src,
             rho_src      = rho_src,
             psim_src     = psim_src,
             psip_src     = psip_src,
             verbosity    = verbosity,
             rho_f =rho_f,
             u_f = u_f,
             E_f = E_f,
             gamma_value = gamma_value,
             cv_value = cv_value,
             check_balance = True)

          # compute exact hydro solution
          hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
             rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

          rad_exact = computeAnalyticRadSolution(mesh, t_end,psim=psim_f,psip=psip_f)
 
          #Compute error
          err.append(computeHydroL2Error(hydro_new, hydro_exact, rad_new, rad_exact ))

          n_elems  *= 2
          dt_value *= 0.5

          # compute convergence rates
          rates_dx = computeHydroConvergenceRates(dx,err)
          rates_dt = computeHydroConvergenceRates(dt,err)

          # print convergence table
          if n_cycles > 1:
             printHydroConvergenceTable(dx,err,rates=rates_dx,
                dx_desc='dx',err_desc='$L_2$')
             printHydroConvergenceTable(dt,err,rates=rates_dt,
                dx_desc='dt',err_desc='$L_2$')

      # plot
      if __name__ == '__main__':

         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)


         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),exact=hydro_exact)

         #plot exact and our E_r
         Er_exact_fn = 1./GC.SPD_OF_LGT*(psim + psip)
         Fr_exact_fn = (psip - psim)*RU.mu["+"]
         Er_exact = []
         Fr_exact = []
         psip_exact = []
         psim_exact = []
         x = mesh.getCellCenters()
         for xi in x:
             
             substitutions = {'x':xi, 't':t_end}
             Er_exact.append(Er_exact_fn.subs(substitutions))
             Fr_exact.append(Fr_exact_fn.subs(substitutions))
             psip_exact.append(psip_f(xi,t_end))
             psim_exact.append(psim_f(xi,t_end))

         plotRadErg(mesh, rad_new.E, Fr_edge=rad_new.F, exact_Er=Er_exact, exact_Fr =
               Fr_exact)

         plotS2Erg(mesh, rad_new.psim, rad_new.psip, exact_psim=psim_exact,
                 exact_psip=psip_exact)

         plotTemperatures(mesh, rad_new.E, hydro_states=hydro_new, print_values=True)

         #Make a pickle to save the error tables
         from sys import argv
         pickname = "results/testRadHydroDiffMMS.pickle"
         if len(argv) > 2:
            if argv[1] == "-o":
               pickname = argv[2].strip()

         #Create dictionary of all the data
         big_dic = {"dx": dx}
         big_dic["dt"] =  dt
         big_dic["Errors"] = err
         pickle.dump( big_dic, open( pickname, "w") )

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
          hydro_BC = HydroBC(bc_type='periodic', mesh=mesh)
      
          # create cross sections
          cross_sects = [(ConstantCrossSection(sig_s, sig_s+sig_a),
                          ConstantCrossSection(sig_s, sig_s+sig_a))
                          for i in xrange(mesh.n_elems)]

          # transient options
          t_start  = 0.0
          t_end = 0.001

          # if run standalone, then be verbose
          if __name__ == '__main__':
             verbosity = 2
          else:
             verbosity = 0

          #slope limiter
          limiter = 'none'
          
          # run the rad-hydro transient
          rad_new, hydro_new = runNonlinearTransient(
             mesh         = mesh,
             problem_type = 'rad_hydro',
             dt_option    = 'constant',
             dt_constant  = dt_value,
             slope_limiter = limiter,
             use_2_cycles = True,
             t_start      = t_start,
             t_end        = t_end,
             rad_BC       = rad_BC,
             cross_sects  = cross_sects,
             rad_IC       = rad_IC,
             hydro_IC     = hydro_IC,
             hydro_BC     = hydro_BC,
             mom_src      = mom_src,
             E_src        = E_src,
             rho_src      = rho_src,
             psim_src     = psim_src,
             psip_src     = psip_src,
             verbosity    = verbosity,
             rho_f =rho_f,
             u_f = u_f,
             E_f = E_f,
             gamma_value = gamma_value,
             cv_value = cv_value,
             check_balance = True)

          # compute exact hydro solution
          hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

          #Compute error
          err.append(computeHydroL2Error(hydro_new, hydro_exact))

          dt_value /= 2.

      # compute convergence rates
      rates = computeHydroConvergenceRates(dt,err)

      # plot
      if __name__ == '__main__':

         # plot radiation solution

         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

         # print convergence table
         if n_cycles > 1:
            printHydroConvergenceTable(dt,err,rates=rates,
               dx_desc='dt',err_desc='L2')

         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),
            exact=hydro_exact)

         #plot exact and our E_r
         Er_exact_fn = 1./GC.SPD_OF_LGT*(psim + psip)
         Fr_exact_fn = (psip - psim)*RU.mu["+"]
         Er_exact = []
         Fr_exact = []
         psip_exact = []
         psim_exact = []
         x = mesh.getCellCenters()
         for xi in x:
             
             substitutions = {'x':xi, 't':t_end}
             Er_exact.append(Er_exact_fn.subs(substitutions))
             Fr_exact.append(Fr_exact_fn.subs(substitutions))
             psip_exact.append(psip_f(xi,t_end))
             psim_exact.append(psim_f(xi,t_end))


         plotRadErg(mesh, rad_new.E, Fr_edge=rad_new.F, exact_Er=Er_exact, exact_Fr =
               Fr_exact)

         plotRadErg(mesh, rad_new.psim, rad_new.psip, exact_Er=psip_exact,
                 exact_Fr=psim_exact)

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
      psim = psim.subs(substitutions)
      psip = psip.subs(substitutions)
      rho_f  = lambdify((symbols('x'),symbols('t')), rho,  "numpy")
      u_f    = lambdify((symbols('x'),symbols('t')), u,    "numpy")
      mom_f  = lambdify((symbols('x'),symbols('t')), mom,  "numpy")
      E_f    = lambdify((symbols('x'),symbols('t')), E,    "numpy")
      psim_f = lambdify((symbols('x'),symbols('t')), psim, "numpy")
      psip_f = lambdify((symbols('x'),symbols('t')), psip, "numpy")
      
      # create uniform mesh
      n_elems = 50
      width = 1.0
      mesh = Mesh(n_elems, width)

      # compute radiation IC
      psi_IC = computeRadiationVector(psim_f, psip_f, mesh, t=0.0)
      rad_IC = Radiation(psi_IC)

      # create rad BC object
      rad_BC = RadBC(mesh, 'dirichlet',psip_BC=psip_f,psim_BC=psim_f)

      # compute hydro IC
      hydro_IC = computeAnalyticHydroSolution(mesh,t=0.0,
         rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

      # create hydro BC
      hydro_BC = HydroBC(bc_type='periodic', mesh=mesh)
  
      # create cross sections
      cross_sects = [(ConstantCrossSection(sig_s, sig_s+sig_a),
                      ConstantCrossSection(sig_s, sig_s+sig_a))
                      for i in xrange(mesh.n_elems)]

      # transient options
      t_start  = 0.0
#      t_end = 0.005
      t_end = 0.02

      # if run standalone, then be verbose
      if __name__ == '__main__':
         verbosity = 2
      else:
         verbosity = 0

      #slope limiter
      limiter = 'none'
      
      # run the rad-hydro transient
      rad_new, hydro_new = runNonlinearTransient(
         mesh         = mesh,
         problem_type = 'rad_hydro',
         dt_option    = 'CFL',
         CFL          = 0.5,
       #  dt_option    = 'constant',
       #  dt_constant  = 0.0002,
         slope_limiter = limiter,
         time_stepper = 'BDF2',
         use_2_cycles = True,
         t_start      = t_start,
         t_end        = t_end,
         rad_BC       = rad_BC,
         cross_sects  = cross_sects,
         rad_IC       = rad_IC,
         hydro_IC     = hydro_IC,
         hydro_BC     = hydro_BC,
         mom_src      = mom_src,
         E_src        = E_src,
         rho_src      = rho_src,
         psim_src     = psim_src,
         psip_src     = psip_src,
         verbosity    = verbosity,
            rho_f =rho_f,
            u_f = u_f,
            E_f = E_f,
            gamma_value = gamma_value,
            cv_value = cv_value,
         check_balance = True)

      # plot
      if __name__ == '__main__':

         # plot radiation solution

         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),exact=hydro_exact)

         #plot exact and our E_r
         Er_exact_fn = 1./GC.SPD_OF_LGT*(psim + psip)
         Er_exact = []
         x = mesh.getCellCenters()
         for xi in x:
             
             substitutions = {'x':xi, 't':t_end}
             Er_exact.append(Er_exact_fn.subs(substitutions))

         plotRadErg(mesh, rad_new.E, exact_Er=Er_exact)

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
      width = 1.0

      # initialize lists for mesh size and L1 error for each cycle
      max_dx = list()
      err = list()

      # loop over refinement cycles
      for cycle in xrange(n_cycles):

         if __name__ == '__main__':
            print("\nCycle %d of %d: n_elems = %d" % (cycle+1,n_cycles,n_elems))

         # create uniform mesh
         mesh = Mesh(n_elems, width)

         # append max dx for this cycle to list
         max_dx.append(mesh.max_dx)
   
         # compute radiation IC
         psi_IC = computeRadiationVector(psim_f, psip_f, mesh, t=0.0)
         rad_IC = Radiation(psi_IC)
   
         #Make rad BC object with vacuum for hydro only
         rad_BC = RadBC(mesh, "vacuum")

         # compute hydro IC
         hydro_IC = computeAnalyticHydroSolution(mesh, t=0.0,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)
   
         # create hydro BC
         hydro_BC = HydroBC(bc_type='dirichlet', mesh=mesh, rho_BC=rho_f,
            mom_BC=mom_f, erg_BC=E_f)

         # create cross sections
         cross_sects = [(ConstantCrossSection(sig_s, sig_s+sig_a),
                         ConstantCrossSection(sig_s, sig_s+sig_a))
                         for i in xrange(mesh.n_elems)]
   
         # if run standalone, then be verbose
         if __name__ == '__main__':
            if n_cycles == 1:
               verbosity = 2
            else:
               verbosity = 1
         else:
            verbosity = 0

         # run the rad-hydro transient
         rad_new, hydro_new = runNonlinearTransient(
            mesh         = mesh,
            problem_type = 'rad_hydro',
            dt_option    = 'CFL',
            CFL          = 0.5,
            slope_limiter = slope_limiter,
            time_stepper = 'BDF2',
            use_2_cycles = True,
            t_start      = 0.0,
            t_end        = t_end,
            rad_BC       = rad_BC,
            hydro_BC     = hydro_BC,
            cross_sects  = cross_sects,
            rad_IC       = rad_IC,
            hydro_IC     = hydro_IC,
            mom_src      = mom_src,
            E_src        = E_src,
            psim_src     = psim_src,
            psip_src     = psip_src,
            rho_src      = rho_src,
            verbosity    = verbosity,
            check_balance = False,
            rho_f        = rho_f,
            u_f          = u_f,
            E_f          = E_f,
            gamma_value  = gamma_value,
            cv_value     = cv_value  )
   
         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)
   
         # compute error
         err.append(computeHydroL2Error(hydro_new, hydro_exact))

         # double number of elements for next cycle
         n_elems *= 2

      # compute convergence rates
      rates = computeHydroConvergenceRates(max_dx,err)

      # print convergence table and plot
      if __name__ == '__main__':

         # print convergence table
         if n_cycles > 1:
            printHydroConvergenceTable(max_dx,err,rates=rates,
               dx_desc='dx',err_desc='L2')

         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),
            exact=hydro_exact)

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
          
          # run the rad-hydro transient
          rad_new, hydro_new = runNonlinearTransient(
             mesh         = mesh,
             problem_type = 'rad_hydro',
             dt_option    = 'constant',
             dt_constant  = dt_value,
             slope_limiter = limiter,
             time_stepper = 'BDF2',
             use_2_cycles = True,
             t_start      = t_start,
             t_end        = t_end,
             rad_BC       = rad_BC,
             cross_sects  = cross_sects,
             rad_IC       = rad_IC,
             hydro_IC     = hydro_IC,
             hydro_BC     = hydro_BC,
             mom_src      = mom_src,
             E_src        = E_src,
             rho_src      = rho_src,
             psim_src     = psim_src,
             psip_src     = psip_src,
             verbosity    = verbosity,
             rho_f        = rho_f,
             u_f          = u_f,
             E_f          = E_f,
             gamma_value  = gamma_value,
             cv_value     = cv_value,
             check_balance = False)

          # compute exact hydro solution
          hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
             rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

          rad_exact = computeAnalyticRadSolution(mesh, t_end,psim=psim_f,psip=psip_f)
 
          #Compute error
          err.append(computeHydroL2Error(hydro_new, hydro_exact, rad_new, rad_exact ))

          n_elems  *= 2
          dt_value *= 0.5

          # compute convergence rates
          rates_dx = computeHydroConvergenceRates(dx,err)
          rates_dt = computeHydroConvergenceRates(dt,err)

          # print convergence table
          if n_cycles > 1:
             printHydroConvergenceTable(dx,err,rates=rates_dx,
                dx_desc='dx',err_desc='$L_2$')
             printHydroConvergenceTable(dt,err,rates=rates_dt,
                dx_desc='dt',err_desc='$L_2$')

      # plot
      if __name__ == '__main__':

         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)


         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),exact=hydro_exact)

         #plot exact and our E_r
         Er_exact_fn = 1./GC.SPD_OF_LGT*(psim + psip)
         Fr_exact_fn = (psip - psim)*RU.mu["+"]
         Er_exact = []
         Fr_exact = []
         psip_exact = []
         psim_exact = []
         x = mesh.getCellCenters()
         for xi in x:
             
             substitutions = {'x':xi, 't':t_end}
             Er_exact.append(Er_exact_fn.subs(substitutions))
             Fr_exact.append(Fr_exact_fn.subs(substitutions))
             psip_exact.append(psip_f(xi,t_end))
             psim_exact.append(psim_f(xi,t_end))

         plotRadErg(mesh, rad_new.E, Fr_edge=rad_new.F, exact_Er=Er_exact, exact_Fr =
               Fr_exact)

         plotS2Erg(mesh, rad_new.psim, rad_new.psip, exact_psim=psim_exact,
                 exact_psip=psip_exact)

         plotTemperatures(mesh, rad_new.E, hydro_states=hydro_new, print_values=True)

         #Make a pickle to save the error tables
         from sys import argv
         pickname = "results/testRadHydroStreamingMMSC100.pickle"
         if len(argv) > 2:
            if argv[1] == "-o":
               pickname = argv[2].strip()

         #Create dictionary of all the data
         big_dic = {"dx": dx}
         big_dic["dt"] =  dt
         big_dic["Errors"] = err
         pickle.dump( big_dic, open( pickname, "w") )

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
      psi_right = psim_f(x=width, t=0.0)

      # initialize lists for mesh size and L1 error for each cycle
      max_dx = list()
      err = list()

      # loop over refinement cycles
      for cycle in xrange(n_cycles):

         if __name__ == '__main__':
            print("\nCycle %d of %d: n_elems = %d" % (cycle+1,n_cycles,n_elems))

         # create uniform mesh
         mesh = Mesh(n_elems, width)

         # append max dx for this cycle to list
         max_dx.append(mesh.max_dx)
   
         # compute radiation IC
         psi_IC = computeRadiationVector(psim_f, psip_f, mesh, t=0.0)
         rad_IC = Radiation(psi_IC)
   
         # compute hydro IC
         hydro_IC = computeAnalyticHydroSolution(mesh, t=0.0,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)
   
         # create hydro BC
         hydro_BC = HydroBC(bc_type='periodic', mesh=mesh)

         # create cross sections
         cross_sects = [(ConstantCrossSection(sig_s, sig_s+sig_a),
                         ConstantCrossSection(sig_s, sig_s+sig_a))
                         for i in xrange(mesh.n_elems)]
   
         # slope limiter option
         slope_limiter = "none"
   
         # if run standalone, then be verbose
         if __name__ == '__main__':
            if n_cycles == 1:
               verbosity = 2
            else:
               verbosity = 1
         else:
            verbosity = 0

         # run the rad-hydro transient
         rad_new, hydro_new = runNonlinearTransient(
            mesh         = mesh,
            problem_type = 'rad_hydro',
            dt_option    = 'CFL',
            CFL          = 0.5,
            slope_limiter = slope_limiter,
            time_stepper = 'BDF2',
            use_2_cycles = True,
            t_start      = t_start,
            t_end        = t_end,
            psi_left     = psi_left,
            psi_right    = psi_right,
            hydro_BC     = hydro_BC,
            cross_sects  = cross_sects,
            rad_IC       = rad_IC,
            hydro_IC     = hydro_IC,
            mom_src      = mom_src,
            E_src        = E_src,
            psim_src     = psim_src,
            psip_src     = psip_src,
            rho_src      = rho_src,
            verbosity    = verbosity,
            check_balance = True,
            rho_f =rho_f,
            u_f = u_f,
            E_f = E_f,
            gamma_value = gamma_value,
            cv_value = cv_value  )
   
         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)
   
         # compute error
         err.append(computeHydroL2Error(hydro_new, hydro_exact))

         # double number of elements for next cycle
         n_elems *= 2

      # compute convergence rates
      rates = computeHydroConvergenceRates(max_dx,err)

      # print convergence table and plot
      if __name__ == '__main__':

         # print convergence table
         if n_cycles > 1:
            printHydroConvergenceTable(max_dx,err,rates=rates,
               dx_desc='dx',err_desc='L2')

         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),
            exact=hydro_exact)

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getCellCenters [as 别名]

#.........这里部分代码省略.........
         display_equations = False)

      # create functions for exact solutions
      substitutions = dict()
      substitutions['alpha'] = alpha_value
      substitutions['c']     = GC.SPD_OF_LGT
      rho = rho.subs(substitutions)
      u   = u.subs(substitutions)
      mom = rho*u
      E   = E.subs(substitutions)
      psim = psim.subs(substitutions)
      psip = psip.subs(substitutions)
      rho_f  = lambdify((symbols('x'),symbols('t')), rho,  "numpy")
      u_f    = lambdify((symbols('x'),symbols('t')), u,    "numpy")
      mom_f  = lambdify((symbols('x'),symbols('t')), mom,  "numpy")
      E_f    = lambdify((symbols('x'),symbols('t')), E,    "numpy")
      psim_f = lambdify((symbols('x'),symbols('t')), psim, "numpy")
      psip_f = lambdify((symbols('x'),symbols('t')), psip, "numpy")
      
      # create uniform mesh
      width = 1.0
      mesh = Mesh(n_elems, width)

      # compute radiation IC
      psi_IC = computeRadiationVector(psim_f, psip_f, mesh, t=0.0)
      rad_IC = Radiation(psi_IC)

      # compute radiation BC; assumes BC is independent of time
      psi_left  = psip_f(x=0.0,   t=0.0)
      psi_right = psim_f(x=width, t=0.0)

      #Create Radiation BC object
      rad_BC = RadBC(mesh, "dirichlet", psi_left=psi_left, psi_right=psi_right)

      # compute hydro IC
      hydro_IC = computeAnalyticHydroSolution(mesh,t=0.0,
         rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

      # create hydro BC
      hydro_BC = HydroBC(bc_type='dirichlet', mesh=mesh, rho_BC=rho_f,
         mom_BC=mom_f, erg_BC=E_f)
  
      # create cross sections
      cross_sects = [(ConstantCrossSection(sig_s, sig_s+sig_a),
                      ConstantCrossSection(sig_s, sig_s+sig_a))
                      for i in xrange(mesh.n_elems)]

      # if run standalone, then be verbose
      if __name__ == '__main__':
         verbosity = 2
      else:
         verbosity = 0

      # run the rad-hydro transient
      rad_new, hydro_new = runNonlinearTransient(
         mesh         = mesh,
         problem_type = 'rad_hydro',
         dt_option    = 'CFL',
         #dt_option    = 'constant',
         CFL          = 0.5,
         #dt_constant  = 0.002,
         slope_limiter = slope_limiter,
         time_stepper = 'BDF2',
         use_2_cycles = True,
         t_start      = 0.0,
         t_end        = t_end,
         rad_BC       = rad_BC,
         cross_sects  = cross_sects,
         rad_IC       = rad_IC,
         hydro_IC     = hydro_IC,
         hydro_BC     = hydro_BC,
         mom_src      = mom_src,
         E_src        = E_src,
         rho_src      = rho_src,
         psim_src     = psim_src,
         psip_src     = psip_src,
         verbosity    = verbosity,
         check_balance = False)

      # plot
      if __name__ == '__main__':

         # compute exact hydro solution
         hydro_exact = computeAnalyticHydroSolution(mesh, t=t_end,
            rho=rho_f, u=u_f, E=E_f, cv=cv_value, gamma=gamma_value)

         # plot hydro solution
         plotHydroSolutions(mesh, hydro_new, x_exact=mesh.getCellCenters(),
            exact=hydro_exact)

         # compute exact radiation energy
         Er_exact_fn = 1./GC.SPD_OF_LGT*(psim + psip)
         Er_exact = []
         x = mesh.getCellCenters()
         for xi in x:
             substitutions = {'x':xi, 't':t_end}
             Er_exact.append(Er_exact_fn.subs(substitutions))

         # plot radiation energy
         plotRadErg(mesh, rad_new.E, exact_Er=Er_exact)