Python Mesh.getElement方法代码示例

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getElement [as 别名]
   def test_computeL1Error1D(self):
      # create mesh
      left = 2.5
      width = 1.5
      mesh = Mesh(3, width, x_start=left)

      # create "numerical solution" data
      numerical_solution = [(1.0,2.0),(1.6,2.2),(2.5,2.7)]

      # specify "exact solution" function
      def exact(x):
         return x**2 + 3.0

      # compute exact integral of difference
      exact_integral = 0.0
      for i in xrange(mesh.n_elems):
         # express local numerical solution as linear function y(x) = m*x + b
         el = mesh.getElement(i)
         xL = el.xl
         xR = el.xr
         yL = numerical_solution[i][0]
         yR = numerical_solution[i][1]
         dx = xR - xL
         dy = yR - yL
         m = dy / dx
         b = yL - xL*m

         # compute local integral of difference
         local_integral = (xR**3 - xL**3)/3.0 - 0.5*m*(xR**2 - xL**2)\
            + (3.0-b)*(xR - xL)
 
         # add to global integral of difference
         exact_integral += local_integral

      # compute numerical integral of difference
      numerical_integral = computeL1ErrorLD(mesh, numerical_solution, exact)

      # assert that numerical and exact integrals are approximately equal
      n_decimal_places = 14
      self.assertAlmostEqual(numerical_integral,exact_integral,n_decimal_places)

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getElement [as 别名]
   def test_TRTBE(self):

      # create mesh
      n_elems = 200
      mesh = Mesh(n_elems,1.)

      # time step size and transient start and end times
      dt      = 0.001
      t_start = 0.0
      #t_end   = 0.01
      t_end   = 0.1

      # initialize temperature
      T_init = 0.05
      T_l    = 0.5
      T_r    = 0.05

      # gamma constant
      gam = 1.4
 
      # material 1 properties and IC
      sig_s1 = 0.0
      sig_a1 = 0.2
      c_v1   = convSpecHeatErgsEvToJksKev(1.E+12)
      rho1   = 0.01
      e1     = T_init*c_v1

      # material 2 properties and IC
      sig_s2 = 0.0
      sig_a2 = 2000.
      c_v2   = c_v1
      rho2   = 10.
      e2     = T_init*c_v2

      # construct cross sections and hydro IC
      cross_sects = list()
      hydro_IC = list()
      for i in range(mesh.n_elems):  
      
         if mesh.getElement(i).x_cent < 0.5: # material 1
            cross_sects.append( (ConstantCrossSection(sig_s1, sig_s1+sig_a1),
                                 ConstantCrossSection(sig_s1, sig_s1+sig_a1)) )
            hydro_IC.append(
               HydroState(u=0,rho=rho1,e=e1,spec_heat=c_v1,gamma=gam))
         else: # material 2
            cross_sects.append((ConstantCrossSection(sig_s2, sig_a2+sig_s2),
                                ConstantCrossSection(sig_s2, sig_a2+sig_s2)))
            hydro_IC.append(
               HydroState(u=0,rho=rho2,e=e2,spec_heat=c_v2,gamma=gam))
  
      # create hydro BC
      hydro_BC = HydroBC(bc_type='reflective', mesh=mesh)
  
      # initialize radiation to equilibrium solution
      psi_left  = computeEquivIntensity(T_l)
      psi_right = computeEquivIntensity(T_r)
      rad_IC    = Radiation([psi_right for i in range(n_elems*4)])
      rad_BC    = RadBC(mesh, "dirichlet", psi_left=psi_left, psi_right=psi_right)

      # time-stepper
      time_stepper = "BDF2"

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

      # run transient
      rad_new, hydro_new = runNonlinearTransient(
         mesh         = mesh,
         time_stepper = time_stepper,
         problem_type = 'rad_mat',
         dt_option    = 'constant',
         dt_constant  = dt,
         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,
         verbosity    = verbosity,
         check_balance = True)

      # plot solutions if run standalone
      if __name__ == "__main__":
          plotTemperatures(mesh, rad_new.E, hydro_states=hydro_new, print_values=False)

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

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

          # mesh
          width = 2.*pi
          mesh = Mesh(n_elems,width)

          # compute hydro IC for the sake of computing initial time step size

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

          # compute the initial time step size according to CFL conditon (actually
          # half). We will then decrease this time step by the same factor as DX each
          # cycle
          if cycle == 0:

              sound_speed = [sqrt(i.p * i.gamma / i.rho) + abs(i.u) for i in hydro_IC]
              dt_vals = [cfl_value*(mesh.elements[i].dx)/sound_speed[i]
                 for i in xrange(len(hydro_IC))]
              dt_value = min(min(dt_vals),0.5) #don't take too big of time step
            
              print "initial dt_value", dt_value

              #Adjust the end time to be an exact increment of dt_values
              print "old t_end: ", t_end
              print "new t_end: ", t_end

          print "This cycle's dt value: ", dt_value

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

          dt.append(dt_value)
          dx.append(mesh.getElement(0).dx)

          # 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, 'periodic')

          # 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)

          # Dimensionless parameters. These are all evaluated at peaks of trig
          # functions, very hard coded
          print "---------------------------------------------"
          print " Diffusion limit info:"
          print "---------------------------------------------"
          print "Size in mfp of cell", mesh.getElement(0).dx*sig_t
          print "Ratio of radiation energy to kinetic", Er(0.75,0)/(rho_f(0.25,0)*u_f(0.0,0)**2), GC.RAD_CONSTANT*T_inf**4/(rho_inf*a_inf**2)
          print "Ratio of speed of light to material sound speed", GC.SPD_OF_LGT/u_f(0.0,0), GC.SPD_OF_LGT/(a_inf)
          print "---------------------------------------------"

          # 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

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getElement [as 别名]
def solveHydroProblem():

    # option to print solution
    print_solution = False

    #-------------------------------------------------------------------------------
    # Construct initial hydro states, probably create an initializer class
    #-------------------------------------------------------------------------------
    width = 1.0
    t_end = 0.05

    t = 0.0
    cfl = 0.5
    n = 500

    #Left and right BC and initial values
    gamma = 1.4 #gas constant
    p_left = 1.0
    u_left  = .750
    rho_left = 1.

    p_right = 0.1
    u_right = 0.0
    rho_right = 0.125

    #made up spec heat
    spec_heat = 1.0

    #Create a mesh, currently hardcoded
    mesh = Mesh(n, width)
    dx = mesh.getElement(0).dx

    if n % 2 != 0: #simple error raise, % has lots of diff meanings in python
        raise ValueError("Must be even number of cells")

    i_left = int(0.3*n)
    i_right = int(0.7*n)

    #Create cell centered variables
    states_a = [HydroState(u=u_left,p=p_left,gamma=gamma,rho=rho_left,
       spec_heat=spec_heat) for i in range(i_left)]
    states_a = states_a + [HydroState(u=u_right,p=p_right,gamma=gamma,
       spec_heat=spec_heat,rho=rho_right) for i in range(i_right)]

    # initialize BC object
    bc = HydroBC(bc_type='reflective', mesh=mesh)

    #-----------------------------------------------------------------
    # Solve Problem
    #----------------------------------------------------------------

    #loop over time steps
    time_index = 0
    while (t < t_end):

        # increment time index
        time_index += 1

        #Compute a new time step size based on CFL
        c = [sqrt(i.p*i.gamma/i.rho)+abs(i.u) for i in states_a]
        dt_vals = [cfl*(mesh.elements[i].dx)/c[i] for i in range(len(c))]
        dt = min(dt_vals)

        t += dt
        #shorten last step to be exact
        if (t > t_end):
            t -= dt
            dt = t_end - t + 0.000000001
            t += dt

        print("Time step %d: t = %f -> %f" % (time_index,t-dt,t))

        # update boundary values
        bc.update(states=states_a, t=t-dt)

        # compute slopes
        slopes = HydroSlopes(states_a, bc=bc)

        #Solve predictor step
        states_a = hydroPredictor(mesh, states_a, slopes, dt)

        # update boundary values
        bc.update(states=states_a, t=t-0.5*dt)

        #Solve corrector step
        states_a = hydroCorrector(mesh, states_a, dt, bc=bc)


    # plot solution
    plotHydroSolutions(mesh,states=states_a)

    # print solution
    if print_solution:
       for state in states_a:
          print state

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

#.........这里部分代码省略.........
          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_value, sig_s_value+sig_a_value),
                          ConstantCrossSection(sig_s_value, sig_s_value+sig_a_value))
                          for i in xrange(mesh.n_elems)]

          # compute the initial time step size according to CFL conditon (actually
          # half). We will then decrease this time step by the same factor as DX each
          # cycle
          if cycle == 0:

              sound_speed = [sqrt(i.p * i.gamma / i.rho) + abs(i.u) for i in hydro_IC]
              dt_vals = [cfl_value*(mesh.elements[i].dx)/sound_speed[i]
                 for i in xrange(len(hydro_IC))]
              dt_value = min(dt_vals)
              print "dt_value for hydro: ", dt_value

              #Make sure not taking too large of step for radiation time scale
              period = 2.*math.pi/Ctilde
              dt_value = min(dt_value, period/8.)
            
              print "initial dt_value", dt_value

              #Adjust the end time to be an exact increment of dt_values
              print "t_end: ", t_end

          print "This cycle's dt value: ", dt_value
          dt.append(dt_value)
          dx.append(mesh.getElement(0).dx)

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

          #slope limiter
          limiter = 'double-minmod'
          
          # 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,

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

#.........这里部分代码省略.........
         print "Temperature", T1, T2
         print "momentum", rho1*u1, rho2*u2
         print "E",E1, E2
         print "E_r",Erad_left, Erad_right

         # final time
         t_end = 1.0

         # temperature plot filename
         test_filename = "radshock_marcho3.pdf"

         # temperature plot exact solution filename
         exact_solution_filename = None


      else:
         raise NotImplementedError("Invalid test case")
         
      # compute radiation BC; assumes BC is independent of time
      c = GC.SPD_OF_LGT
      psi_left  = 0.5*c*Erad_left  
      psi_right = 0.5*c*Erad_right

      #Create BC object
      rad_BC = RadBC(mesh, "dirichlet", psi_left=psi_left, psi_right=psi_right)
      
      # construct cross sections and hydro IC
      cross_sects = list()
      hydro_IC = list()
      psi_IC = list()

      for i in range(mesh.n_elems):  
      
         if mesh.getElement(i).x_cent < mesh_center: # material 1
            cross_sects.append( (ConstantCrossSection(sig_s1, sig_s1+sig_a1),
                                 ConstantCrossSection(sig_s1, sig_s1+sig_a1)) )
            hydro_IC.append(
               HydroState(u=u1,rho=rho1,e=e1,spec_heat=c_v1,gamma=gam))

            psi_IC += [psi_left for dof in range(4)]

         else: # material 2
            cross_sects.append((ConstantCrossSection(sig_s2, sig_a2+sig_s2),
                                ConstantCrossSection(sig_s2, sig_a2+sig_s2)))
            hydro_IC.append(
               HydroState(u=u2,rho=rho2,e=e2,spec_heat=c_v2,gamma=gam))

            psi_IC += [psi_right for dof in range(4)]

      #Smooth out the middle solution optionally, shouldnt need this
      n_smoothed = 0
      state_l = hydro_IC[0]
      state_r = hydro_IC[-1]
        
      rho_l =  state_l.rho
      rho_r =  state_r.rho
      drho  = rho_r - rho_l
      u_l   =  state_l.u
      u_r   =  state_r.u
      du    = u_r - u_l
      e_l   =  state_l.e
      e_r   =  state_r.e
      de    = e_r-e_l
      print "p", state_l.p, state_r.p
      print "mach number", state_l.u/state_l.getSoundSpeed(), state_r.u/state_r.getSoundSpeed()

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getElement [as 别名]
   def test_Hydro(self):

      # create mesh
      n_elems = 100
      width = 1.0
      mesh = Mesh(n_elems,width)
      x_diaphragm = 0.3

      # time step size and transient start and end times
      CFL     = 0.5
      t_start = 0.0
      t_end   = 0.2

      # slope limiter
      slope_limiter = "vanleer"

      # constant properties
      sig_s = 1.0 # arbitrary
      sig_a = 0.0 # set to zero to ensure no emission
      sig_t = sig_s + sig_a
      c_v   = 3.0 # arbitrary
      gam = 1.4

      # hydro IC values for left half of domain
      rhoL = 1.0
      uL   = 0.75
      pL   = 1.0

      # hydro IC values for right half of domain
      rhoR = 0.125
      uR   = 0.0
      pR   = 0.1

      # compute left and right internal energies
      eL = pL/(rhoL*(gam - 1.0))
      eR = pR/(rhoR*(gam - 1.0))

      # construct cross sections and hydro IC
      cross_sects = list()
      hydro_IC = list()
      for i in range(mesh.n_elems):  

         # cross section is constant
         cross_sects.append( (ConstantCrossSection(sig_s, sig_t),
                              ConstantCrossSection(sig_s, sig_t)) )
      
         # IC for left half of domain
         if mesh.getElement(i).x_cent < x_diaphragm:
            hydro_IC.append(
               HydroState(u=uL,rho=rhoL,e=eL,gamma=gam,spec_heat=c_v) )

         # IC for right half of domain
         else:
            hydro_IC.append(
               HydroState(u=uR,rho=rhoR,e=eR,gamma=gam,spec_heat=c_v) )

      # create hydro BC
      hydro_BC = HydroBC(bc_type='reflective', mesh=mesh)
  
      # initialize radiation to zero solution to give pure hydrodynamics
      rad_IC    = Radiation([0.0 for i in range(n_elems*4)])
      psi_left  = 0.0
      psi_right = 0.0
      rad_BC    = RadBC(mesh, "dirichlet", psi_left=psi_left, psi_right=psi_right)

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

      # run transient
      rad_new, hydro_new = runNonlinearTransient(
         mesh         = mesh,
         time_stepper = 'BE',
         problem_type = 'rad_hydro',
         dt_option    = 'CFL',
         CFL          = 0.5,
         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,
         slope_limiter = slope_limiter,
         verbosity    = verbosity,
         check_balance= True)



      # plot solutions if run standalone
      if __name__ == "__main__":

         # get exact solutions
         #get the exact values
         f = open('exact_testHydro.txt', 'r')
         x_e = []
         u_e = []
#.........这里部分代码省略.........

# 需要导入模块: from mesh import Mesh [as 别名]
# 或者: from mesh.Mesh import getElement [as 别名]
   def test_RadHydroShock(self):

      # create uniform mesh
      n_elems = 500
      width = 0.04
      x_start = -0.02
      mesh_center = x_start + 0.5*width
      mesh = Mesh(n_elems, width, x_start=x_start)

      # slope limiter
      slope_limiter = "double-minmod"

      # gamma constant
      gam = 5.0/3.0
 
      # material 1 and 2 properties;
      T_ref = 0.1 #keV reference unshocked, upstream, ambient, equilibrium
      
      sig_a1 = 577.35
      sig_s1 = 0.0
      c_v1   = 0.14472799784454
      sig_a2 = sig_a1
      sig_s2 = sig_s1
      c_v2   = c_v1

      #Read in Jim's nondimensional results to set preshock and postshock dimensional
      mach_number = "1.2" #Choices are 2.0, 1.2, 3.0, 5.0
      filename = 'analytic_shock_solutions/data_for_M%s.pickle' % mach_number
      f = open(filename,'r')
      data = pickle.load(f)
      f.close()

      #compute scalings based on an assumed density and reference temperature
      dp =getDimensParams(T_ref=T_ref, rho_ref=1.0, C_v=c_v1, gamma=gam)
      
      print data

      #Scale non-dimensional values into dimensional results
      rho1   = data['Density'][0]*dp['rho']
      u1     = data['Speed'][0]*dp['a']   #velocity times reference sound speed
      Erad1  = data['Er'][0]*dp['Er']
      T1     = data['Tm'][0]*dp["Tm"]
      e1     = T1*c_v1
      E1     = rho1*(e1 + 0.5*u1*u1)

      # material 2 IC
      rho2   = data['Density'][-1]*dp['rho']
      u2     = data['Speed'][-1]*dp['a']   #velocity times reference sound speed
      Erad2  = data['Er'][-1]*dp['Er']
      T2     = data['Tm'][-1]*T_ref
      e2     = T2*c_v2
      E2     = rho2*(e2 + 0.5*u2*u2)
 
      print r"$\rho$ & %0.8e & %0.8e & g cm$^{-3}$ \\" % (rho1, rho2)
      print r"$u$ & %0.8e & %0.8e & cm sh$^{-1}$ \\" % (u1, u2)
      print r"$T$ & %0.8e & %0.8e & keV \\" % (T1,T2)
      print r"$E$ & %0.8e & %0.8e & Jks cm$^{-3}$\\" % (E1,E2)
      print r"$E_r$ & %0.8e & %0.8e & Jks cm$^{-3}$ \\" % (Erad1, Erad2)
      print r"$F_r$ & %0.8e & %0.8e & Jks cm$^{-2}$ s$^{-1}$ \\" %(0.0,0.0)
      print r"vel", u1, "&", u2
      print "Temperature", T1,"&",  T2
      print "momentum", rho1*u1,"&",  rho2*u2
      print "E",E1,"&",  E2
      print "E_r",Erad1, "&", Erad2
      print sig_a1, sig_s1
      exit()
      # material 1 IC
 
      # final time
      t_end = 0.8
 
      # temperature plot filename
      test_filename = "radshock_mach_"+re.search("M(\d\.\d)",filename).group(1)+".pdf"
 
      # compute radiation BC; assumes BC is independent of time
      c = GC.SPD_OF_LGT
      psi_left  = 0.5*c*Erad1  
      psi_right = 0.5*c*Erad2

      #Create BC object
      rad_BC = RadBC(mesh, "dirichlet", psi_left=psi_left, psi_right=psi_right)
      
      # construct cross sections and hydro IC
      cross_sects = list()
      hydro_IC = list()
      psi_IC = list()

      for i in range(mesh.n_elems):  
      
         if mesh.getElement(i).x_cent < mesh_center: # material 1
            cross_sects.append( (ConstantCrossSection(sig_s1, sig_s1+sig_a1),
                                 ConstantCrossSection(sig_s1, sig_s1+sig_a1)) )
            hydro_IC.append(
               HydroState(u=u1,rho=rho1,e=e1,spec_heat=c_v1,gamma=gam))

            psi_IC += [psi_left for dof in range(4)]

         else: # material 2
            cross_sects.append((ConstantCrossSection(sig_s2, sig_a2+sig_s2),
                                ConstantCrossSection(sig_s2, sig_a2+sig_s2)))
#.........这里部分代码省略.........