35.7 Other Test Problems

35.7.1 The non-equilibrium ionization test problem

The neitest problem tests the ability of FLASH to calculate non-equilibrium ionization (NEI) ion abundances. It simulates a stationary plasma flow through a temperature step profile. The solutions were checked using an independent stationary code based on a fifth order Runge-Kutta method with adaptive stepsize control by step-doubling (see Orlando et al. (1999)).

Figure 35.74: Temperature profile assumed for the test.

The test assumes a plasma with a mass density of $2\times 10^{-16}$ gm cm$^{-3}$ flowing with a constant uniform velocity of $3\times 10^5$ cm s$^{-1}$ through a temperature step between $10^4$ K and $10^6$ K (cf. Figure 35.74). The plasma is in ionization equilibrium before going through the jump in the region at $T = 10^4$ K. The population fractions in equilibrium are obtained from the equations

$$[n_{i}^Z]_{eq} S_{i}^Z = [n_{i + 1}^Z]_{eq} \alpha_{i + 1}^Z  (i = 1, ..., l_Z - 1)$$ (35.55)

$$\sum_{i=1}^{l_Z} [n_{i}^Z]_{eq} = A_{Z} n_{p}$$ (35.56)

The presence of a temperature jump causes a strong pressure difference, which in turn should cause significant plasma motions. Since the purpose is to test the NEI module, it is imposed that the pressure difference does not induce any plasma motion and, to this end, the hydro variables (namely, $T$, $\rho$, v) are not updated. In practice, the continuity equations are solved with uniform density and velocity, while the momentum and energy equations are ignored.

Figure 35.75 shows the population fractions for the 12 most abundant elements in astrophysical plasmas derived with the stationary code (Orlando et al. (1999)). The out of equilibrium ionization conditions are evident for all the elements just after the flow goes through the temperature jump.

Figure 35.75: Numerical solutions of the stationary code. The figure shows the population fractions vs. space for the 12 elements most abundant in astrophysical plasmas assuming a stationary flow through a temperature jump.

Figure 35.75: ... continued ...

Figure: As in Figure 35.75 for the solutions of the FLASH code.

Figure 35.76: ... continued ...

The same problem was solved with the NEI module of the FLASH code, assuming that the plasma is initially in ionization equilibrium at $t=t_0$ over all the spatial domain. After a transient lasting approximately 700 s, in which the population fractions evolve due to the plasma flow through the temperature jump, the system reaches the stationary configuration. Outflow boundary conditions (zero-gradient) are assumed at both the left and right boundaries. Figure 35.76 shows the population fraction vs. space after 700 s.

35.7.2 The Delta-Function Heat Conduction Problem

The ConductionDelta problem tests the basic function of the Diffuse unit, in particular its use for thermal conduction, in connection with the Conductivity material property unit. It can be easily modified to examine the effects of Viscosity, of non-constant conductivity, and of combining these diffusive effects with hydrodynamics.

In its default configuration, ConductionDelta models the time development of a temperature peak that at some time $t_{\rm ini}$ has the shape of a delta function in 1D or 3D, subject to heat conduction (with constant coefficient) only. An ideal gas EOS is assumed. When using a flux-based Diffuse interface, the setup includes the Hydro code unit, but changes of the solution due to hydrodynamical effects are completely suppressed by zeroing all hydrodynamic fluxes each time before diffusive fluxes are computed (using $\htmladdnormallink{{\tt updateHydroFluxes}}{http://flash.rochester.edu/site/fla... ...txt}\index[rp]{Hydro!updateHydroFluxes@{\tt updateHydroFluxes}} = {\tt .FALSE.}$); diffusive fluxes are then computed either by calling Diffuse_therm (in the +splitHydro case) or by an calling an internal routine (such as hy_uhd_addThermalFluxes) of the Hydro unit.

The theoretical solution of this initial value problem is known: For any $t > t_{\rm ini}$, the temperature profile is a Gaussian that spreads self-similarly over time. In particular in 1D, if the initial condition is defined as

$$T(x, t_{\rm ini}) = Q\delta(x) \quad,$$ (35.57)

then

$$T(x, t) = \frac{Q}{(4\pi\chi (t-t_{\rm ini}))^{1/2}} e^{-x^2/4\chi (t-t_{\rm ini})}$$ (35.58)

(with $\chi$ the constant coefficient of diffusivity), see for example Zel'dovich and Raizer Ch. X.

Figure: Temperature profile of the Delta-Function Heat Conduction Problem at two times. $t_{\rm ini}=-1 {\rm ms}$, top: $t=0 {\rm ms}$, bottom: $t=2.46 {\rm ms}$.

See the end of Sec:Diffuse for alternative ways of configuring the test problem, using either a flux-based or a standalone Diffuse interface.

35.7.3 The HydroStatic Test Problem

The Hydrostatic problem tests the basic function of hydrostatic boundary conditions implemented in the Grid unit, in connection with a Hydro implementation. It is essentially a 1D problem, but can be configured as 1 1D, 2D, or 3D setup. It can be easily modified to include additional physical effects, by including additional code units in the setup.

In its default configuration, HydroStatic is set up with constant Gravity. The domain is initialized with density, pressure, etc., fields representing an analytical solution of the hydrostatic problem with the given gravitational acceleration, essentially following the barometric formula.

This initial condition is then evolved in time. Ideally, the solution would remain completely static, and nothing should move. The deviation from this ideal behavior that occurs in practice is a measure of the quality of the discretization of the initial condition, of the hydrodynamics implementation, and of the boundary conditions. The effect of the latter, in particular, can be examined by visualizing artifacts that develop near the boundaries (in particular, velocity artifacts), and studying their dependence on the choice of boundary condition variant.

35.7.4 Hybrid-PIC Test Problems

A classic plasma model problem is that of an ion beam into a plasma (Winske 1986). Work by Matthew (1994) describes a two-dimensional simulation of a low density ion beam through a background plasma. The initial condition has a uniform magnetic field, with a beam of ions, number density $n_{b}$, propagating along the field with velocity $v_b=10v_A$, where $v_A=B/\sqrt{\mu_0nm_i}$ is the Alfvén velocity, through a background (core) plasma of number density $n_{c}$. Both the background and beam ions have thermal velocities $v_{th}=v_A$. Here we study what is denoted a resonant beam with $n_{b}=0.015 n_c$ and $v_c=-0.2 v_A$. Electron temperature is assumed to be zero, so the electron pressure term in the electric field equation of state is zero. The weight of the macroparticles are chosen such that there is an equal number of core and beam macroparticles, each beam macroparticles thus represent fewer real ions than the core macroparticles. The number of magnetic field update subcycles is nine.

35.7.4.1 One-dimensional ion beam

The spatial extent of the domain is 22016 km along the $x$-axis, divided up in 256 cells with periodic boundary conditions. The core number density is $n_c=7$ cm$^{-3}$. The magnetic field magnitude is $B=6$ nT directed along the $x$-axis, which give an Alfvén velocity of 50 km/s and an ion inertial length of $\delta_i$=86 km, where $\delta_i=v_A/\Omega_i$ and $\Omega_i=q_iB/m_i$ is the ion gyrofrequency. The time step is 0.0865 s = $0.05 \Omega_i^{-1}$, and the cell size is $\delta_i$. The number of particles per cell is 32. In Fig. 35.78 we show a velocity space plot of the macro-ions at time $t=34.6$ s $\approx 20 \Omega_i^{-1}$. This can be compared to Fig. 5 in Winske (1984).

Figure: Velocity space plot for a one-dimensional ion beam. Velocity along the $y$-axis as a function of position at time $t=34.6$ s $\approx 20 \Omega_i^{-1}$. Each gray dot is a core macro-ion, and each black dot is a beam macro-ion.

35.7.4.2 Two-dimensional ion beam

In the two-dimensional case we have a square grid with sides of length 22016 km, and 128 cells in each direction with periodic boundary conditions. The time step is 0.0216 s = $0.05 \Omega_i^{-1}$, and the cell widths are $2\delta_i$. The number of particles per cell is 16. Otherwise the setup is identical to the one-dimensional case. In Fig. 35.79 we show the magnitude of the magnetic field $y$-component at time $t=77.85$ s $\approx 46 \Omega_i^{-1}$. This can be compared to Fig. 5 in Winske (1986).

Figure: The magnetic field $y$-component for a two-dimensional ion beam at time $t=77.85$ s $\approx 46 \Omega_i^{-1}$.

35.7.5 Full-physics Laser Driven Simulation

The LaserSlab simulation provides an example of how to run a realistic, laser driven, simulation of an HEDP experiment. It exercises a wide range of HEDP physics including:

• 3T unsplit hydrodynamics in $R-z$ geometry
• Electron thermal conduction with Spitzer conductivities
• Ion/electron equilibration
• Radiation Diffusion with tabulated opacity
• Laser Ray Tracing
• Tabulated Equation of State

This simulation has no analytic solutions associated with it. Rather, it is designed to demonstrate the HEDP capabilities that are available in FLASH.

The simulation models a single laser beam illuminating a solid Aluminum disc in $R-z$ geometry. The laser is focused on the $z$-axis and enters the domain at a 45 degree angle. The simulation contains two materials (or “species” in FLASH terminology). These are called CHAM (short for chamber) and TARG (short for target). The TARG species models the Aluminum disc. The laser travels through a low density Helium region before reaching the Aluminum. The properties of the CHAM species are set to model the low density Helium. A schematic showing the initial conditions is shown in Figure 35.80. The laser rays enter just above the lower-right corner of the domain and travel until they are absorbed or reflected. The beam intensity is not uniform across the beam cross section. Rather, it has a super-Gaussian intensity profile with a 10 micron spot radius.

Figure 35.80: Schematic showing the initial conditions in the LasSlab simulation.

The following setup line is used for this simulation:

  ./setup -auto LaserSlab -2d +cylindrical -nxb=16 -nyb=16 +hdf5typeio \
          species=cham,targ +mtmmmt +laser +uhd3t +mgd mgd_meshgroups=6 \
          -parfile=example.par

where:

• -2d and +cylindrical tells FLASH that this is a 2D cylindrical simulation
• -nxb=16 and -nyb=16 sets the AMR block size to 16 by 16
• +hdf5typeio enables a newer version of IO that writes to HDF5 files. HDF5TypeIO is faster that other IO implementations in FLASH and is also more user friendly when generating figures using external tools, such as VisIt.
• species=cham,targ Creates two species named CHAM and TARG using the setup argument instead of the SPECIES keyword in the Simulation Config file. This is required for using the Multispecies opacity unit and is more user friendly since it allows many options to be specified in the runtime parameters file. See Sec:DefineSpecies2 for more information.
• +mtmmmt Activates the Multi-Temperature, Multi-Material, Multi-Type EOS in FLASH. This EOS lets users select a different EOS model (such as gamma-law or tabulated) for each species in the simulation through the runtime parameters file. (see Sec:3TEos Multitype)
• +laser Activates laser ray tracing (see Sec:EnergyDeposition).
• +mgd and mgd_meshgroups=6 activates multigroup radiation diffusion and tells FLASH that there will be maximum of six groups per mesh (see chp:RadTrans)
• +uhd3t Activates the 3T unsplit hydro solver (see Sec:3THydro)
• -parfile=example.par Chooses a parameter file that is different than the default

The remainder of this section will describe how each physics unit functions and how the runtime parameters are chosen for this simulation. Images of the results will also be shown at the end of this section.

35.7.5.1 Initial Conditions and Material Properties

The section “Initial Conditions” section in the runtime parameter file, “example.par”, sets parameters which define the initial state of the fluid and also set the properties of the materials. The first runtime parameters in this section are:  

sim_targetRadius = 200.0e-04
sim_targetHeight = 20.0e-04
sim_vacuumHeight = 60.0e-04

These options define the space occupied by the target material (the targ species, Aluminum for this simulation). The Aluminum disc extends from:

$$0 < R < {\tt sim\_targetRadius}\;\;\mathrm{and}$$ (35.59)

$${\tt sim\_vacuumHeight} < z < {\tt sim\_vacuumHeight} + {\tt sim\_targetHeight},$$ (35.60)

where all parameters are in centimeters.

The next section sets the initial fluid state in the target region and also sets the properties of the targ species to match Aluminum. The relevant options are:  

# Target material defaults set for Aluminum at room temperature:
sim_rhoTarg  = 2.7                   # Set target material density (solid Al)
sim_teleTarg = 290.11375             # Set initial electron temperature
sim_tionTarg = 290.11375             # Set initial ion temperature 
sim_tradTarg = 290.11375             # Set initial radiation temperature
ms_targA = 26.9815386                # Set average atomic mass
ms_targZ = 13.0                      # Set average atomic number
ms_targZMin = 0.02                   # Set minimum allowed ionization level
eos_targEosType = "eos_tab"          # Set EOS model to tabulated EOS
eos_targSubType = "ionmix4"          # Set EOS table file to IONMIX 4 format
eos_targTableFile = "al-imx-003.cn4" # Set tabulated EOS file name

The initial density is chosen to match solid Aluminum. The initial ion, electron, and radiation temperatures are set to room temperature (290 K, 0.025 eV). The runtime parameters ms_targA and ms_targZ set the average atomic mass and atomic number of the targ species. Note that in the past, this information was specified by modifying the Simulation_initSpecies subroutine. This information can now be specified using runtime parameters since the species=cham,targ setup variable was used (see Sec:DefineSpecies2 for more information). The parameter ms_targZMin sets the minimum ionization level for targ. Without this parameter, the initial ionization would be nearly zero and the slab would be subcritical with respect to the laser light. The laser rays would penetrate deeply into the slab and could possible travel through it with no heating. This can be avoided by putting a lower limit on the ionization level. The final three parameters control the behavior of the MTMMMT EOS with respect to targ. The option eos_targEosType sets the EOS model to use for the Aluminum. In this case, we have chosen to use a tabulated EOS which is in the IONMIX4 format. Finally, we specify the name of the file which contains the EOS data. The file “al-imx-003.cn4” can be found in the LaserSlab simulation directory in the source tree. The IONMIX4 format, described in Sec:OpacityIonmix, is not human readable. However, the file “al-imx-003.imx.gz” contains a human readable representation of this EOS and opacity data.

A similar set of runtime parameter options exist for the cham species:  

# Chamber material defaults set for Helium at pressure 1.6 mbar:
sim_rhoCham  = 1.0e-05
sim_teleCham = 290.11375
sim_tionCham = 290.11375
sim_tradCham = 290.11375
ms_chamA = 4.002602
ms_chamZ = 2.0
eos_chamEosType = "eos_tab"
eos_chamSubType = "ionmix4"
eos_chamTableFile = "he-imx-005.cn4"

The initial Helium density, sim_rhoCham is chosen to be fairly low. This ensures that the Helium remains relatively transparent to the laser light (see Sec:LaserPowerDeposition for a description of how the laser energy absorption coefficient is computed).

For both the targ and cham species, the initial state of the fluid is defined using a density with three temperatures. The initial EOS mode is set using the eosModeInit runtime parameter. It is set to “dens_temp_gather”. This tells FLASH that the initial state will be specified using the density and the temperature.

35.7.5.2 Multigroup Radiation Diffusion Parameters

This simulation uses multigroup radiation diffusion (MGD) to model the radiation field. This involves setting parameters which control the MGD package itself and setting parameters which control how opacities for the two species are computed. The MGD settings from the parameter file are:  

rt_useMGD       = .true.
rt_mgdNumGroups = 6
rt_mgdBounds_1  = 1.0e-01
rt_mgdBounds_2  = 1.0e+00
rt_mgdBounds_3  = 1.0e+01
rt_mgdBounds_4  = 1.0e+02
rt_mgdBounds_5  = 1.0e+03
rt_mgdBounds_6  = 1.0e+04
rt_mgdBounds_7  = 1.0e+05
rt_mgdFlMode    = "fl_harmonic"
rt_mgdFlCoef    = 1.0


rt_mgdXlBoundaryType = "reflecting"
rt_mgdXrBoundaryType = "vacuum"
rt_mgdYlBoundaryType = "vacuum"
rt_mgdYrBoundaryType = "reflecting"
rt_mgdZlBoundaryType = "reflecting"
rt_mgdZrBoundaryType = "reflecting"

The first option turns on MGD. The second parameter sets the number of radiation energy groups in the simulation, six in this case. Note that this number must be less than or equal to the value of the mgd_meshgroups setup variable. This restriction can be relaxed through the use of mesh replication (see Sec:MgdMeshRep). The runtime parameters beginning with rt_mgdBounds_ define the group boundaries for each energy group in electron-volts. The parameters rt_mgdFlMode and rt_mgdFlCoef set the flux limiter coefficient (see Sec:DiffuseFluxLimiters and Sec:MGD for more information). The final six parameters set the boundary conditions for all groups.

The next set of parameters tell FLASH how to compute opacities for each material. In general, the total cell opacity is a number density weighted average of the opacity of each species in the cell. The parameters which control the behavior of the opacity unit for this simulation are:  

useOpacity     = .true.


### SET CHAMBER (HELIUM) OPACITY OPTIONS ###
op_chamAbsorb   = "op_tabpa"
op_chamEmiss    = "op_tabpe"
op_chamTrans    = "op_tabro"
op_chamFileType = "ionmix4"
op_chamFileName = "he-imx-005.cn4"


### SET TARGET (ALUMINUM) OPACITY OPTIONS ###
op_targAbsorb   = "op_tabpa"
op_targEmiss    = "op_tabpe"
op_targTrans    = "op_tabro"
op_targFileType = "ionmix4"
op_targFileName = "al-imx-003.cn4" 

The first parameter tells FLASH that opacities will be computed. The next five parameters control how opacities are computed for the cham species (Helium) and the final five control the targ species (Aluminum). These options tell FLASH to use the IONMIX4 formatted opacity file “he-imx-005.cn4” for Helium and “al-imx-003.cn4” for Aluminum. Sec:OpacityMultispecies provides a detailed description of how these parameters are used. Note that the number of groups and group structure in these opacity files must be consistent with the rt_mgdNumGroups and rt_mgdBounds parameters that were described above.

35.7.5.3 Laser Parameters

This section describes the runtime parameters which control the behavior of the laser in FLASH. Modeling laser energy deposition using ray tracing is fairly complicated and requires large numbers of input parameters. It is highly recommended that you read through the section describing the behavior of the laser model in FLASH to learn more about how the ray tracing algorithms actually work (see Sec:EnergyDeposition).

This simulation uses a single laser beam which travels at an angle of 45 degrees from the $z$-axis. The “Laser Parameters” section of the runtime parameters file contains all of the options relevant to defining the laser beam. The first set of parameters are:  

useEnergyDeposition = .true. # Turn on laser energy deposition
ed_maxRayCount      = 10000  # Set maximum number of rays/cycle/proc to 2000
ed_gradOrder        = 2      # Turn on laser ray refraction

The first parameter activates the laser model.

The ed_maxRayCount parameter sets the maximum number of rays which can be stored on a given process at once. It tells FLASH how much space to set aside for storing rays. In the case of a single process simulation, ed_maxRayCount must be less than the total number of rays launched on each cycle. The total number of rays launched is set by the ed_numRays_### parameters, described below. Finally, the parameter ed_gradOrder controls how the electron number density and temperature are interpolated within a cell. A value of two specifies linear interpolation.

The next set of parameters are:  

ed_laser3Din2D           = .true.
ed_laser3Din2DwedgeAngle = 0.1

Setting ed_laser3Din2D to .true. activates the 3D-in-2D ray trace algorithm. This means that in this simulation, the beams are defined in 3D geometry and the laser rays are traced in 3D. The total energy deposited is then projected back on the $R-z$ plane. This can significantly improve the accuracy of the ray trace when compared to 2D-in-2D ray-tracing. The exception to this is the case of a single beam centered on the $z$-axis and traveling along the $z$ axis. In this case, there is no difference between 3D-in-2D and 2D-in-2D ray tracing algorithms. In this example, the rays travel at a 45 degree angle. Thus, the 3D-in-2D ray-trace is needed. When using the 3D-in-2D ray trace, ed_laser3Din2DwedgeAngle sets the wedge angle. It should be set to a number between 0.1 and 1.0 degrees. See Sec:3Din2D.

The next set of parameters are:  

ed_useLaserIO             = .true.
ed_laserIOMaxNumPositions = 10000
ed_laserIONumRays         = 128

These parameters make use of the LaserIO package in FLASH for visualizing laser rays (see Sec:LaserIO). These options tell FLASH to write out the trajectory of 128 rays to the plot files. The runtime parameter plotFileIntervalStep is set to 100 which tells FLASH to write plot files every 10 cycles. Later in this section, instructions and examples for plotting the laser ray data will be presented.

Next, a single laser pulse is defined. Each laser pulse (in this case, there is only one) represents a specific time history of the power and multiple pulses can be defined. The LaserSlab simulation has only a single laser beam which uses laser pulse number 1. The relevant runtime parameter options are:  

# Define Pulse 1:
ed_numberOfSections_1 = 4
ed_time_1_1  = 0.0
ed_time_1_2  = 0.1e-09
ed_time_1_3  = 1.0e-09
ed_time_1_4  = 1.1e-09


ed_power_1_1 = 0.0
ed_power_1_2 = 1.0e+09
ed_power_1_3 = 1.0e+09
ed_power_1_4 = 0.0

The laser pulse is defined as a piecewise linear function. The parameter ed_numSections_1 is set to four, meaning the laser power is defined using four different time/power points. The four time and power pairs are defined using the ed_time_1_# and ed_power_1_# parameters where the times are in seconds and the powers are specified in Watts. In this case a square pulse is defined with a gradual rise, termination over 0.1 ns. The peak power is $10^9$ W and the pulse lasts for a total of 1.1 ns.

The next set of parameters define the laser beam itself. The trailing _1 in the parameter names indicate that these parameters are associated with the first beam (in this case, the only beam):  

ed_lensX_1                    =  1000.0e-04
ed_lensY_1                    =  0.0e-04
ed_lensZ_1                    = -1000.0e-04
ed_lensSemiAxisMajor_1        =  10.0e-04
ed_targetX_1                  =  0.0e-04
ed_targetY_1                  =  0.0e-04
ed_targetZ_1                  =  60.0e-04
ed_targetSemiAxisMajor_1      =  10.0e-04
ed_targetSemiAxisMinor_1      =  10.0e-04
ed_pulseNumber_1              =  1
ed_wavelength_1               =  1.053
ed_crossSectionFunctionType_1 = "gaussian2D"
ed_gaussianExponent_1         =  4.0
ed_gaussianRadiusMajor_1      =  7.5e-04
ed_gaussianRadiusMinor_1      =  7.5e-04
ed_numberOfRays_1             =  4096
ed_gridType_1                 = "radial2D"
ed_gridnRadialTics_1          =  64
ed_semiAxisMajorTorsionAngle_1=  0.0
ed_semiAxisMajorTorsionAxis_1 = "x"

The parameters ed_targetX_1, ed_targetY_1, and ed_targetZ_1 define the coordinate of the center of the laser focal spot. Notice that even though this is a 2D simulation, three coordinates are specified. That is because this while most of the FLASH solvers will operate in 2D, $R-z$ geometry, the laser ray trace operates in 3D Cartesian geometry. For specifying the laser focal spot center, “X” refers to the $R$ direction in the FLASH simulation (which is also called the “X” direction for many FLASH runtime parameters, regardless of the geometry). The ed_targetZ_1 parameter is referring to the $z$ direction in the $R-z$ simulation. This is actually the “Y” direction for other FLASH runtime parameters. Finally, the “Z” direction for the 3D-in-2D ray trace parameters is the direction pointing out of the computer screen. This same naming convention is used for the parameters ed_lensX_1, ed_lensY_1, and ed_lensZ_1. These parameters specify the center of the lens. On each time step, all rays are traced from a location on the lens to the focal spot. The lens must be defined so that it exists entirely outside of the domain. The parameters ed_targetSemiAxisMajor_1 and ed_targetSemiAxisMinor_1 set the radius of the laser spot. The laser beam can have an elliptical shape. In this case the two parameters will have different values, but usually, a circular beam is desired. In this case, the beam will have a circular shape and a 10 micron radius. The parameter ed_lensSemiAxisMajor_1 sets the radius of the focal spot. The laser beam is associated with pulse number 1 (the only pulse in this case) using the parameter ed_pulseNumber_1. The parameter ed_wavelength_1 specifies the wavelength, in microns, of the laser light. The ed_numberOfRays_1 parameter tells FLASH to launch 4096 rays for the first beam every time the EnergyDeposition subroutine is called. When using unsplit hydrodynamics (as in this case), this routine is called a single time per cycle. When using split hydrodynamics, it is called twice.

Taken together, the parameters:

• ed_crossSectionFunctionType_1,
• ed_gaussianExponent_1,
• ed_gaussianRadiusMajor_1, and
• ed_gaussianRadiusMinor_1

define the spatial variation of the intensity across the beam. When ed_crossSectionFunctionType_1 is set to “gaussian2D”, the intensity profile is a supergaussian described by:

$$I(r) = I_0 \exp \left\{ - \left[ \left( \frac{x}{R_x} \right)^2 + \left( \frac{y}{R_y} \right)^2 \right]^{\gamma} \right\}$$ (35.61)

where:

• $R_x =$ed_gaussianRadiusMajor_1,
• $R_y =$ed_gaussianRadiusMinor_1, and
• $\gamma =$ed_gaussianExponent_1.

In our case, we've set up a beam with a super-Gaussian profile ($gamma = 4$) and an e-folding length of 7.5 microns. Since $R_x = R_y$ we have a circular beam and not an elliptical beam. The user does not directly specify $I_0$. Rather, the user specifies the total beam power (as described above). This then sets $I_0$ through the integral relation:

$$P = 2 \pi \int_0^R dr r I(r),$$ (35.62)

The parameters ed_semiAxisMajorTorsionAngle_1 and ed_semiAxisMajorTorsionAxis_1 only need to change for non-circular beams.

Finally, the user must specify how rays will be spatially distributed across the beam. The ed_numberOfRays_1 parameter sets the total number of rays to launch on each time step for each beam. Now we have to decide how to distribute those rays cross the cross-section of the laser beam. The parameter ed_gridType_1 is set to “radial2D”. This means the rays will be laid out on the circular cross-section of the beam with some radial and angular spacing. The area of the beam is divided into regions of equal radius and angle. The parameter ed_gridnRadialTics_1 is set to 64 meaning that there are 64 radial slices in the beam cross section. Therefore, there are: $4096/64 = 64$ angular slices.

35.7.5.4 Examining LaserSlab Simulation Output

This section describes how the output can be visualized using VisIt and through other means. LaserIO code in FLASH can be used to visualize ray paths (see Sec:LaserIO for more information). The LaserIO unit can only be used with parallel HDF5 IO implementations. The +hdf5typeio setup option is consistent with this. The runtime parameters:  

ed_useLaserIO             = .true.
ed_laserIOMaxNumPositions = 10000
ed_laserIONumRays         = 128

tell FLASH to use LaserIO and to write out the trajectory of 128 rays to the plot files so they can be visualized in VisIt. Notice that this is a small subset of the total number of rays that are actually launched. Writing out all of the rays is inadvisable since it will have numerous negative effects on performance when large number of rays are used in the simulation. In this case, plot files are written every 0.1 ns and the simulation is run for 2 ns cycles (as indicated by the parameter tmax). Thus, the following plot files are generated:  

lasslab_hdf5_plt_cnt_0000
lasslab_hdf5_plt_cnt_0001
...
lasslab_hdf5_plt_cnt_0199
lasslab_hdf5_plt_cnt_0200

These HDF5 files contain the ray trajectories for each of the plot time steps (cycles 1, 11, ... 191, and 201 in this case). Unfortunately, at this time, the FLASH VisIt plugin does not natively support ray visualization. The ray data must be extracted from the plot files in placed in a form that VisIt does understand. The extract_rays.py script, found in the tools/scripts directory, performs this operation. Note, however, that extract_rays.py requires the NumPy and PyTables python packages to operate.

After this example LaserSlab simulation is run, the extract_rays.py script can be used as shown:  

»> extract_rays.py lasslab_hdf5_plt_cnt_*
Processing file: lasslab_hdf5_plt_cnt_0000
Processing file: lasslab_hdf5_plt_cnt_0001
...
Processing file: lasslab_hdf5_plt_cnt_0198
Processing file: lasslab_hdf5_plt_cnt_0199
Processing file: lasslab_hdf5_plt_cnt_0200
No ray data in "lasslab_hdf5_plt_cnt_0200" skipping...

Note that a warning message was printed for the last plot file. This occurs because there is no ray data for a plot file written on the last cycle of a simulation. This is a limitation of the LaserIO package. The script generates files in the VTK format named:  

lasslab_las_0000.vtk
  lasslab_las_0001.vtk
  ...
  lasslab_las_0198.vtk
  lasslab_las_0199.vtk

These files contain meshes which define the ray trajectories in a format understood by VisIt.

The ray data can be plotted on top of the computational mesh. In VisIt, load the lasslab_las_*.vtk database. Add the “RayPower_Watts” Pseudocolor plot. The “RayPower_Watts” plot colors the trajectories based on the power of each ray in units of Watts. Thus, the power of any one ray at any given time should be less than the beam powers specified in the runtime parameters file.

An example showing the ray powers and trajectories is shown in Figure 35.80. The trajectory of 128 laser rays has been drawn. The rays are colored based on their powers.

Computing the Number Density and Average Ionization: For a simulation using the MTMMMT EOS, the electron number density itself is not directly represented in the solution vector. Rather, for mostly historical reasons, FLASH only keeps track of the variables YE and SUMY, which are defined as:

$${\tt YE} = \frac{\bar{Z}}{\bar{A}}\;,\;\; {\tt SUMY} = \frac{1}{\bar{A}}$$

where $\bar Z$ is the average ionization level and $\bar A$ is the average atomic mass of the cell. Thus, one can compute the average ionization and electron number density from these quantities. For example:

$$\bar{A} = \frac{1}{{\tt SUMY}},$$

$$\bar{Z} = \frac{{\tt YE}}{{\tt SUMY}},$$

$$n_\mathrm{ion} = N_A \frac{\rho}{\bar A} = {\tt SUMY} N_A \rho,$$

$$n_\mathrm{ele} = N_A \bar Z \frac{\rho}{\bar A} = {\tt YE} N_A \rho$$

Where $N_A$ is Avogadros number, $n_\mathrm{ion}$ is the ion number density, $n_\mathrm{ele}$ is the electron number density, and $\rho$ is the mass density (DENS in FLASH). These quantities can be plotted in VisIt using expressions.

Another useful tool in visualizing the behavior of the laser is the FLASH variable DEPO. This variable stores the laser energy deposited in a cell on a particular time step per unit mass. The division by mass removes the geometric factor present in the ray power. The amount of laser energy entering the domain and the amount of laser energy exiting the domain can be used to compute the fraction of laser energy that is deposited. The energy entering and exiting the domain is shown for each cycle and integrated over the entire simulation in the file “lasslab_LaserEnergyProfile.dat”. This file is produced whenever the laser ray tracing package is active in FLASH.

35.7.6 Laser Driven Simulation with Thomson Scattering Diagnostics

The LaserSlab simulation described previously can be configured to include the ThomsonScattering diagnostic unit. We describe this configuration here to serve as a demonstration for using that unit. See chp:ThomsonScattering for more information on the diagnostic code.

As in the previously described LaserSlab setup, The simulation models a single laser beam illuminating a solid Aluminum disc in $R-z$ geometry. The laser is focused on the $z$-axis and enters the domain at a 45 degree angle. The simulation contains two materials (or “species” in FLASH terminology). These are called CHAM (short for chamber) and TARG (short for target). The TARG species models the Aluminum disc. The laser travels through a low density Helium region before reaching the Aluminum. The properties of the CHAM species are set to model the low density Helium. A schematic showing the initial conditions is shown in Figure 35.80. The laser rays enter just above the lower-right corner of the domain and travel until they are absorbed or reflected. The beam intensity is not uniform across the beam cross section. Rather, it has a super-Gaussian intensity profile with a 10 micron spot radius.

The following setup line is used for this simulation:

  ./setup -auto LaserSlab -2d -nxb=16 -nyb=16 +hdf5typeio \
          species=cham,targ +mtmmmt +laser +uhd3t +mgd mgd_meshgroups=6 \
          ed_maxBeams=2 ed_maxPulseSections=4 ed_maxPulses=2 \
          thsc_maxBeams=1 thsc_maxPulseSections=4 thsc_maxPulses=1 thsc_maxDetectors=1 \
          -parfile=example_thsc.par ThscDemo=True

Here the setup variable ThscDemo is used to turn on the optional ThomsonScattering configuration of the LaserSlab setup. This variable is specific to the LaserSlab simulation; see the end of this simulation's Config file to see how this is done for LaserSlab and use these ot similar directives to enable ThomsonScattering for other setups.

Additionally,

• some setup variables beginning with thsc_ are used to demonstrate how to set some limits for parameters related to ThomsonScattering that will be hardwired into the executable;
• some setup variables beginning with ed_ are used to demonstrate how to set some limits for parameters related to EnergyDeposition that will be hardwired into the executable;
• -parfile=example_ths.par chooses a parameter file that has additional settings for runtime parameters that control ThomsonScattering;
• other arguments are as for the regular LaserSlab setup and have been described previously.

35.7.7 Z-pinch Simulation

The ZPinch simulation provides an example of how to run a realistic, current-driven, simulation of a Z-pinch experiment. It exercises a wide range of HEDP physics including:

• 3T unsplit magnetohydrodynamics in $R-z$ geometry
• Current drive from the Circuit unit
• Implicit magnetic resistive diffusion with DaviesWen resistivities
• Electron and ion thermal conduction with JiHeld conductivities
• Ion/electron heat exchange
• Radiation Diffusion with tabulated opacity
• Tabulated Equation of State

This simulation has no analytic solutions associated with it. Rather, it is designed to demonstrate some of the HEDP capabilities that are available in FLASH, and serve as a template for users to simulate their own Z-pinch problems.

The following setup line is used for this simulation:

  ./setup -auto magnetoHD/ZPinch -1d +cylindrical +ug +usm3t +mtmmmt \
          species=fill,line,vacu +mgd mgd_meshgroups=6 +hdf5typeio

The included flash.par file contains many runtime parameters for customization. Below is a subset of the parameters used in ZPinch:  

sim_rhoType = 0  #  0: 0 Gaussians


#### DD Fill
sim_fill_dens = 30.e-03
sim_fill_minDens = 1.e-05   # used in rhoType 11 
sim_fill_ctr = 0.           # Gaussian center
sim_fill_sigma = 1.0        # rho  exp(-r^2/(2*sigma^2))
sim_fill_tele = 1160452.21  # 100 eV (preheated)
sim_fill_tion = 1160452.21  # 100 eV
sim_fill_trad = 1160452.21  # 100 eV
sim_fill_maxTemp = 1.e12    # default, effectively no limit
sim_fill_minTemp = 1.e-12   # effectively no limit
ms_fillA = 2.014
ms_fillZ = 1.0
eos_fillEosType = "eos_tab"
eos_fillSubType = "ionmix4"		# set eos table file to ionmix4 format
eos_fillTableFile = "DD-006-imx.cn4"	# set tabulated eos file name


#### Be Liner
sim_line_dens = 1.848
sim_line_minDens = 1.e-06    # also acts as density floor for vacuum
sim_line_ctr = 1.0           # Gaussian center
sim_line_sigma = 0.8         # rho  exp(-r^2/(2*sigma^2))
sim_innerRadius = 0.29       # inner radius of the liner
sim_thickness = 0.058        # thickness of the liner
sim_pert = 0.		      # pseudo-random density perturbation amplitude
sim_seed = 1		      # seed for random number generator (makes sims repeatable)
sim_line_tele = 5802.26105    # 0.5 eV
sim_line_tion = 5802.26105    # 0.5 eV
sim_line_trad = 5802.26105    # 0.5 eV
sim_line_maxTemp = 1.e12      # default, effectively no limit
sim_line_minTemp = 1.e-12     # effectively no limit
ms_lineA = 9.012
ms_lineZ = 4.
eos_lineEosType = "eos_tab"
eos_lineSubType = "ionmix4"		# set eos table file to ionmix4 format
eos_lineTableFile = "Be-006-imx.cn4"	# set tabulated eos file name


#### Vacuum (low density Be)
sim_vacu_dens = 1.e-05         # rho <= this is vacuum
sim_vacu_tele = 5802.26105     # 0.5 eV
sim_vacu_tion = 5802.26105     # 0.5 eV
sim_vacu_trad = 5802.26105     # 0.5 eV
sim_vacu_maxTemp = 5802.26105  # 0.5 eV
sim_vacu_minTemp = 5802.26105  # 0.5 eV
ms_vacuA = 9.012
ms_vacuZ = 4.
eos_vacuEosType = "eos_tab"
eos_vacuSubType = "ionmix4"		# set eos table file to ionmix4 format
eos_vacuTableFile = "Be-006-imx.cn4"	# set tabulated eos file name


#### Current drive
circ_currFile = "current.dat"	    # filename for current input (ns vs. MA)
#circ_voltFile = "Voc_example.dat"  # filename for voltage input (circuit model) (s vs. V)
#circ_cylLength = 1.                # used in 1D sims with circuit model


# circuit model restart parameters
#circ_initIload	= 1.404306536824466661E+07	# these init parameters get used when a restart
#circ_initIstack= 1.828569946357272193E+07	# is detected. User should refer to 'circuit.dat'
#circ_initVc	= 3.834045253898037598E+06	# for the appropriate values to use.
#circ_initFlux	= 1.724316746860757865E-01


#### trajectory output settings
sim_trajOutputInterval = 1.e-10       	# output interval for trajectory.dat
sim_trajOutputIntervalNearStag = 1.e-11	# output interval near stagnation
sim_rNearStag = 0.145              	# fuel radius at which <= defines "near stagnation"

Many of the parameters are self-explanatory, or the comments included in the flash.par file suffice. The parameter sim_rhoType controls which density profile is used for initialization; 0 is the default value and is simply constant densities in each region. Other options are 1 (single Gaussian for fill species), 2 (a double Gaussian for fill and line species, summed), and 11 (a double Gaussian for fill and line species, piecewise). The _ctr and _sigma parameters are used for Gaussian profile options. For more details about how sim_rhoType is used, see the Simulation_initBlock.F90 subroutine.

The _minTemp and _maxTemp parameters provide temperature control for each region. These are enforced every timestep in the Simulation_adjustEvolution.F90 subroutine and are typically most useful for the vacuum region. The Simulation_adjustEvolution.F90 subroutine also outputs a file named trajectory.dat, which contains radial positions and moving-average velocities of the inner and outer liner interfaces. The sim_trajOutputInterval, sim_trajOutputIntervalNearStag, and sim_rNearStag parameters provide control of the output frequency of this trajectory data.

The default Config file is set to use the FileInput implementation of the Circuit unit (see Sec:Circuit. Furthermore, the parameter circ_currFile defines the name of the provided current input file as current.dat. Parameters related to the McBride circuit model implementation are commented out; to use this, the user first needs to change the requested implementation in the Config file, resetup, and recompile.

This simulation uses a high constant resistivity for the vacuum and the DaviesWen resistivity coefficients elsewhere. See vacuum resistivity for details on how to control resistivity and Ohmic heating in a vacuum region.

The default parameters included in FLASH give a 1D Z-pinch that stagnates around 160 ns. The uniform grid settings use 128 cells across the 0.5 cm domain, and iProcs defines the number of processors to run on (currently set to 32). In Figure , radial profiles of the azimuthal magnetic field are shown for three different simulation times. We see that the target pinches due to $\mathbf{J}\times\mathbf{B}$ forces in a reasonable fashion and the vacuum region maintains a $B/r$ profile (zero current).

Figure 35.81: Radial profiles of azimuthal magnetic field for the ZPinch example problem.

There is a 2D version of the flash.par file included, but it likely requires much higher resolution to obtain reasonable results.