diff --git a/docs/source/verification/verification.rst b/docs/source/verification/verification.rst index 42353c9f8..7a317cb8c 100644 --- a/docs/source/verification/verification.rst +++ b/docs/source/verification/verification.rst @@ -9,21 +9,21 @@ decaying Taylor vortex is used. .. math:: :label: advConvTV_u - + u = u_o - cos(\pi(x-u_ot)) sin(\pi(y-v_ot))e^{-2.0\omega t} - + .. math:: :label: advConvTV_v - v = v_o + sin(\pi(x-u_ot)) cos(\pi(y-v_ot))e^{-2.0\omega t} - + v = v_o + sin(\pi(x-u_ot)) cos(\pi(y-v_ot))e^{-2.0\omega t} + .. math:: :label: advConvTV_p - + p = -\frac{p_o}{4}(cos(2\pi(x-u_ot)) + cos(2\pi(y-v_ot)))e^{-4\omega t} - + In this study, the constants :math:`u_o`, :math:`v_o`, and :math:`p_o` are all assigned values of :math:`1.0`, and the viscosity :math:`\mu` is @@ -72,9 +72,9 @@ the full set of advection operators supported in Nalu-Wind. :math:`v` component of velocity using fourth order pressure stabilization with timestep scaling, backward Euler - + .. _so-fourth-tstep: - + .. figure:: figures/convTaylorVortexSO.pdf :width: 500px :align: center @@ -83,9 +83,9 @@ the full set of advection operators supported in Nalu-Wind. :math:`v` component of velocity using fourth order pressure stabilization with timestep scaling, BDF2 - + .. _hybrid-tstep: - + .. figure:: figures/convTaylorVortexSO_ElemLagElemPf.pdf :width: 500px :align: center @@ -178,16 +178,16 @@ manufactured as follows: v &= +v_o sin(a \pi x) cos(a \pi y ) sin(a \pi z) \\ w &= -w_o sin(a \pi x) sin(a \pi y ) cos(a \pi z) \\ p &= -\frac{p_o}{4}( cos(2 a \pi x) + cos(2 a \pi y ) + cos(2 a \pi z) ) \\ - h &= +h_o cos(a_h \pi x) cos(a_h \pi y ) cos(a_h \pi z) + h &= +h_o cos(a_h \pi x) cos(a_h \pi y ) cos(a_h \pi z) The equation of state is simply the ideal gas law, .. math:: :label: ideal-gas-eos1 - + \rho = \frac{P^{ref} M}{R T} - + The simulation is run on a three-dimensional domain ranging from -0.05:+0.05 with constants @@ -219,17 +219,17 @@ manufactured as follows: v &= +v_o sin(a \pi x) cos(a \pi y ) sin(a \pi z) \\ w &= -w_o sin(a \pi x) sin(a \pi y ) cos(a \pi z) \\ p &= -\frac{p_o}{4}( cos(2 a \pi x) + cos(2 a \pi y ) + cos(2 a \pi z) ) \\ - z &= +z_o cos(a_z \pi x) cos(a_z \pi y ) cos(a_z \pi z) - + z &= +z_o cos(a_z \pi x) cos(a_z \pi y ) cos(a_z \pi z) + The equation of state is simply the standard inverse mixture fraction property expression for density, .. math:: :label: ideal-gas-eos2 - + \rho = \frac{1} {\frac{z}{rho^P} + \frac{1-z}{rho^S} } - + The simulation is run on a three-dimensional domain ranging from -0.05:+0.05 with constants :math:`a, a_z, \rho^p, \rho^s, Sc, \mu` equal @@ -270,9 +270,9 @@ temperature field takes on the following form: .. math:: :label: adv-conv-tv-u - + T = \frac{\lambda}{4} (cos(2 a \pi x) + cos(2 a \pi y)). - + The above manufactured solution is run on three meshes of domain size of 1x1. The domain was first meshed as a triangular mesh and then converted @@ -287,9 +287,9 @@ as follows: .. math:: :label: diff-op - + -\int \Gamma \frac{\partial \phi}{\partial x_j} A_j. - + The choice of the gradient operator at the integration point is a functin of the underlying method. For CVFEM, the gradient operator is @@ -297,25 +297,25 @@ provided by the standard shape function derivatives, .. math:: :label: cvfem-derivative2 - + \frac{\partial \phi_{ip}}{\partial x_j} = \sum \frac{\partial N^{ip}_{j,k}} {\partial x_j} \phi_k. - + For the edge-based scheme, a blending of an orthogonal gradient along the edge and a NOC is employed, .. math:: :label: general-grad - \frac{\partial \phi_{ip}}{\partial x_j} = \bar{G_j\phi} + \left[ \left(\phi_R - \phi_L \right) + \frac{\partial \phi_{ip}}{\partial x_j} = \bar{G_j\phi} + \left[ \left(\phi_R - \phi_L \right) - \bar{G_l\phi}dx_l \right] \frac{A_j}{A_k dx_k}. In the above equation, :math:`G_j\phi` is a projected nodal gradient. The general equation for this quantity is .. math:: :label: png - + \int w_I G_j \phi {dV} = \int w_i \frac{\partial \phi}{\partial x_j}{dV}. - + Possible forms of this include either lumped or consistent mass (the later requires a global equation solve) with either the full CVFEM @@ -360,12 +360,12 @@ meshes for the CVFEM operator (R3 and R4). Results indicate that the .. _laplace-tquad-a: -.. figure:: figures/tquadLaplaceMMS.pdf +.. figure:: figures/tquadLaplaceMMS.pdf :width: 500px :align: center - Error norms for tquad4 refinement study. - R0, R1, and R2 refinement. + Error norms for tquad4 refinement study. + R0, R1, and R2 refinement. .. _laplace-tquad-b: @@ -373,7 +373,7 @@ meshes for the CVFEM operator (R3 and R4). Results indicate that the :width: 500px :align: center - Error norms for tquad4 refinement study. + Error norms for tquad4 refinement study. R0, R1, R2, R3, R4, and R4 refinementError for CVFEM. An inspection of the magnitude of error between the exact and computed @@ -399,9 +399,9 @@ evaluate the DG-based CVFEM approach. .. math:: :label: threed-t - + T = \frac{\lambda}{4} (cos(2 a \pi x) + cos(2 a \pi y) + cos(2 a \pi z)). - + Figure :numref:`laplace-nc` represents the MMS field for temperature. The simulation study includes uniform refinement of a first- and @@ -430,7 +430,7 @@ Figure :numref:`laplace-ncoc-b`. MMS order of accuracy for nonconformal algorithm. Temperature norms for P1 and P2 elements. - + .. _laplace-ncoc-b: .. figure:: figures/dgNonconformalCVFEM_3dTempMMS_OoCPNG.pdf @@ -463,13 +463,13 @@ results for a P1 CVFEM implementation. Precursor-based Simulations ------------------------------------------------------ In the field of turbulent flow modeling and simulation, often times simulations -may require sophisticated boundary conditions that can only be obtained +may require sophisticated boundary conditions that can only be obtained from previously run high-fidelity simulations. For example, consider a typical turbulent jet simulation in which the experimental inlet condition was preceeded by a turbulent pipe entrance region. Furthermore, in most cases the ability to adequately predict the developing jet flow regime may be highly sensitive to proper inlet conditions. Figure :numref:`inlet-pipe` and -Figure :numref:`inlet-pipe-jet` outline a process in which a high fidelity +Figure :numref:`inlet-pipe-jet` outline a process in which a high fidelity large-eddy simulation of a periodic pipe was used to determine a representative inlet condition for a turbulent round jet. Specifically, a precursor pipe flow simulation is run with velocity provided to an output file. This output file serves @@ -493,35 +493,35 @@ as the inlet velocity profile for the subsequent simulation. Subsequent turbulent jet simulation using the precursor data obtained by a periodic pipe flow. -In the above use case, as with most general simulation studies, the mesh resolution -for the precursor simulation may be different from the subsequent simulation. Moreover, +In the above use case, as with most general simulation studies, the mesh resolution +for the precursor simulation may be different from the subsequent simulation. Moreover, the time scale for the precursor simulation may be much shorter than the subsequent simulation. Finally, the data required for the subsequent simulation will likely be at different time steps unless an overly restrictive rule is enforced, i.e., a fixed timestep -for each simulation. +for each simulation. -In order to support such use cases, extensive usage of the the Sierra Toolkit infrastructure -is expected, most notably within the IO and Transfer modules. The IO module can be used to interpolate -the precursor simulation boundary data to the appropriate time required by the subsequent +In order to support such use cases, extensive usage of the the Sierra Toolkit infrastructure +is expected, most notably within the IO and Transfer modules. The IO module can be used to interpolate +the precursor simulation boundary data to the appropriate time required by the subsequent simulation. Specifically, the IO module linearly interpolates between the closest data interval in the precursor data set. A recycling offset factor is included within the IO interface that allows -for the cycling of data over the full time scale of interest within the subsequent simulation. For +for the cycling of data over the full time scale of interest within the subsequent simulation. For typical statistically stationary turbulent flows, this is useful to ensure proper statistics are captured in subsequent runs. -After the transient data set from the precursor simulation is interpolated to the proper time, -the data is spatially interpolated and transferred to the subsequent simulation mesh using the -STK Transfer module. Efficient coarse parallel searches (point/bounding box) provide the list of -candidate owning elements on which the fine-scale search operates to determine the best search -candidate. The order of spatial interpolation depends on the activated numerical discretization. +After the transient data set from the precursor simulation is interpolated to the proper time, +the data is spatially interpolated and transferred to the subsequent simulation mesh using the +STK Transfer module. Efficient coarse parallel searches (point/bounding box) provide the list of +candidate owning elements on which the fine-scale search operates to determine the best search +candidate. The order of spatial interpolation depends on the activated numerical discretization. Therefore, by combining the two STK modules, the end use case to support data transfers of boundary data is supported. As noted, there are many other use cases in addition to the overviewed turbulent jet simulation -that require such temporal/spatial interpolation capabilities. For example, in typical wind +that require such temporal/spatial interpolation capabilities. For example, in typical wind farm simulation applications, a proper atmospheric boundary layer (ABL) configuration is required -to capture a given energy state of the boundary layer. In this case, a periodic precusor ABL is run -with the intent of providing the inlet condition to the subsequent wind farm domain. As with the +to capture a given energy state of the boundary layer. In this case, a periodic precusor ABL is run +with the intent of providing the inlet condition to the subsequent wind farm domain. As with the previous description, the infrustructure requirements remain the same. Finally, the general creation of an "input_output" region can be useful in validation cases @@ -530,9 +530,9 @@ PLIF experimental data sets. Although the temporal interpolation is not required of this data at high time step frequency is useful for post-processing. In this verification section, a unit test approach will be referenced that is captured within the -STK module test suite. This foundational test coverage provides confidence in the underlying IO and +STK module test suite. This foundational test coverage provides confidence in the underlying IO and parallel search/interpolation processes. In addition to briefly describing the infrastructure testing, -application tests are provided as further evidence of correctness. The application test first is based +application tests are provided as further evidence of correctness. The application test first is based on the convecting Taylor vortex verification case while the second is the ABL precursor application space demonstration. @@ -540,8 +540,8 @@ Infrastructure Unit Test ++++++++++++++++++++++++ As noted above, the Nalu-Wind application code leverages the STK unit tests within the IO and transfer modules. Interested parties may peruse the STK product under a cloned Trilinos cloned project, -i.e., Trilinos/packages/stk/stk_doc_test. Under the STK product, a variety of search, transfer and -input/output tests exist. For example, interpolation in time using the IO infrastructure is captured +i.e., Trilinos/packages/stk/stk_doc_test. Under the STK product, a variety of search, transfer and +input/output tests exist. For example, interpolation in time using the IO infrastructure is captured in addition to a variety of search and transfer use cases. Application Verification Test; Convecting Taylor Vortex @@ -555,7 +555,7 @@ a very fine mesh simulation is run with boundary conditions specified in the inp to be of type, "convecting_taylor_vortex". This specifies the analytical function for the x-component of velocity as provided in Equation :eq:`advConvTV_u`. The simulation is run while providing output to a Realm of type "input_output" using a transfer objective, -"input_output". The transient data is then used +"input_output". The transient data is then used for a series of mesh refinement studies. The viscosity is set at 0.001 while the domain is 1x1. In this study, the edge-based scheme is activated, however, the precursor interpolation methodology is not sensitive to the underlying numerical method. @@ -580,7 +580,7 @@ are represented as the subsets of the full data - both in space and time. As suc degradation of second-order accuracy is expected. The subsequent simulations are run with an "external_data" transfer objective and a Realm of type, "external_data_provider". -In Figure :numref:`ctv-l2`, a plot of :math:`L_2` norms of the x-component of velocity are shown +In Figure :numref:`ctv-l2`, a plot of :math:`L_2` norms of the x-component of velocity are shown for the subsequent set of simulations that use the precursor data. Results of this study verify both the second-order temporal accuracy of the underlying numerical scheme and the process of using both space and time interpolation. @@ -600,17 +600,17 @@ Application Verification Test; ABL Precursor/Subsequent The second, and final application test is an ABL-based simulation that first runs a precursor periodic solution in order to capture an appropriate ABL specification. The boundary data saved from the precursor -simulation are then used as an inflow boundary condition for the subsequent ABL simulation. As the precurosr -is run for a smaller time frame than the subsequent simulation, the usage of data cycling is active. This -full integration test is captured within the regression test suite. The simulation is described as a non-isothermal +simulation are then used as an inflow boundary condition for the subsequent ABL simulation. As the precurosr +is run for a smaller time frame than the subsequent simulation, the usage of data cycling is active. This +full integration test is captured within the regression test suite. The simulation is described as a non-isothermal turbulent flow. In Figure :numref:`abl-susequent-cycle`, the transient recycling of the ABL thermal inflow boundary -condition is captured at an arbitrary nodal location very near the wall boundary condition. The subsequent -simulation reads the precursor data set for time zero seconds until 3000 seconds at which time it recylces -the inlet condition back to the initial precursor simulation time, i.e., zero seconds. An interesting note in this -study is the fact that the precursor periodic simulation, which was run at the same Courant number, was using time -steps approximately three times greater than the subsequent inflow/open configuration. +condition is captured at an arbitrary nodal location very near the wall boundary condition. The subsequent +simulation reads the precursor data set for time zero seconds until 3000 seconds at which time it recylces +the inlet condition back to the initial precursor simulation time, i.e., zero seconds. An interesting note in this +study is the fact that the precursor periodic simulation, which was run at the same Courant number, was using time +steps approximately three times greater than the subsequent inflow/open configuration. .. _abl-susequent-cycle: @@ -621,9 +621,9 @@ steps approximately three times greater than the subsequent inflow/open configur Transient recycling of the temperature inflow boundary condition for the subsequent ABL simulation. After 3000 seconds, the inflow boundary condition is recycled from the begining of the precursor simulation. -In Figure :numref:`abl-susequent-check-one-two`, (left) the subsequent simulation inflow temperature field and -full profile over the full domain is captured at approximately 4620 seconds. On the right of the figure, the -temperature boundary condition data that originated from the precursor simulation, which was read into the subsequent +In Figure :numref:`abl-susequent-check-one-two`, (left) the subsequent simulation inflow temperature field and +full profile over the full domain is captured at approximately 4620 seconds. On the right of the figure, the +temperature boundary condition data that originated from the precursor simulation, which was read into the subsequent "external_field_provider" Realm, is shown (again at approximately 4620 seconds). .. _abl-susequent-check-one-two: @@ -632,7 +632,7 @@ temperature boundary condition data that originated from the precursor simulatio :width: 500px :align: center - Subsequent simulation showing the full temperature domain (left) and on the precursor inflow temperature boundary + Subsequent simulation showing the full temperature domain (left) and on the precursor inflow temperature boundary condition field obtained from the perspective of the subsequent "external_field_provider" Realm (right). @@ -642,17 +642,15 @@ Boussinesq Verification Unit tests ++++++++++ -Unit-level verification was performed for the Boussinesq body force term :eq:`boussbuoy` with a -nodal source appropriate to the edge-based scheme (MomentumBoussinesqSrcNodeSuppAlg.single_value) as well as a -separate unit test for the element-based "consolidated" Boussinesq source term -(MomentumKernelHex8Mesh.buoyancy_boussinesq). Proper volume integration with different element topologies is +Unit-level verification was performed for the Boussinesq body force term :eq:`boussbuoy` with a +nodal source. Proper volume integration with different element topologies is also tested (the "volume integration" tests in the MasterElement and HOMasterElement test cases). Stratified MMS +++++++++++++++++++++++++++++++++++++++++++++++++++++++ -A convergence study using the method of manufactured solutions (MMS) was also performed to assess the integration +A convergence study using the method of manufactured solutions (MMS) was also performed to assess the integration of the source term into the governing equations. An initial condition of a Taylor-Green vortex for velocity, a zero- gradient pressure field, and a linear enthalpy profile in the z-direction are imposed. @@ -666,14 +664,14 @@ gradient pressure field, and a linear enthalpy profile in the z-direction are im h &= z. The simulation is run on a three-dimensional domain ranging from -1/2:+1/2 with reference density, -reference temperature and the thermal expansion coefficient to equal to 1, 300, and 1, respectively. -:math:`\beta` is much larger than typical (:math:`1 / T_{\rm ref}`) so that the buoyancy term is a +reference temperature and the thermal expansion coefficient to equal to 1, 300, and 1, respectively. +:math:`\beta` is much larger than typical (:math:`1 / T_{\rm ref}`) so that the buoyancy term is a significant term in the MMS in this configuration. -The Boussinesq buoyancy model uses a gravity vector of magnitude of ten in the z-direction +The Boussinesq buoyancy model uses a gravity vector of magnitude of ten in the z-direction opposing the enthalpy gradient, :math:`g_i = (0, 0, -10)^T`. The temperature for this test ranges -between 250K and 350K. The test case was run with a regular hexahedral mesh, using the edge-based -vertex centered finite volume scheme. Each case was run with a fixed maximum Courant number of 0.8 +between 250K and 350K. The test case was run with a regular hexahedral mesh, using the edge-based +vertex centered finite volume scheme. Each case was run with a fixed maximum Courant number of 0.8 relative to the specified solution. @@ -733,7 +731,7 @@ relative to the specified solution. +---------------+---------------------+---------------+---------------+-------+ -This test is added to Nalu-Wind's nightly test suite, testing that the convergence rate between +This test is added to Nalu-Wind's nightly test suite, testing that the convergence rate between the 1/32 and 1/64 element sizes is second order. 3D Hybrid 1x2x10 Duct: Specified Pressure Drop @@ -742,10 +740,10 @@ In this section, a specified pressure drop in a simple 1x2x10 configuration is r a variety of homogeneous blocks of the following topology: hexahedral, tetrahedral, wedge, and thexahedral. This analytical solution is given by an infinite series and is coded as the "1x2x10" user function. The simulation is run with an outer wall boundary condition -with two open boundary conditions. The specified pressure drop is 0.016 over the 10 cm +with two open boundary conditions. The specified pressure drop is 0.016 over the 10 cm duct. The density and viscosity are 1.0e-3 and 1.0e-4, respectively. The siumulation study is run a fixed Courant numbers with a mesh spacing ranging from 0.2 to 0.025. -Figure :numref:`specified-dp-flow-hex-tet` provides the standard velocity profile for the +Figure :numref:`specified-dp-flow-hex-tet` provides the standard velocity profile for the structured hexahedral and unstructured tetrahedral element type. .. _specified-dp-flow-hex-tet: @@ -759,10 +757,10 @@ structured hexahedral and unstructured tetrahedral element type. The simulation study employed a variety of elemental topologies of uniform mesh spacing as noted above. Figure :numref:`specified-dp-l2` outlines the convergence in the :math:`L_2` norm using the low-order elemental CVFEM implementation using the recently changed tetrahedral -and wedge element quadrature rules. Second-order accuracy is noted. Interestingly, the +and wedge element quadrature rules. Second-order accuracy is noted. Interestingly, the hexahedral and wedge topology provided nearly the same accuracy. Also, the tetrahedral accuracy was approximately four tiomes greater. Finally, the Thexahedral topology -proved to be second-order, however, provided very poor accuracy results. +proved to be second-order, however, provided very poor accuracy results. .. _specified-dp-l2: @@ -774,12 +772,12 @@ proved to be second-order, however, provided very poor accuracy results. 3D Hybrid 1x1x1 Cube: Laplace ----------------------------- -The standard Laplace operator is evalued on the full set of low-order hybrid topologies +The standard Laplace operator is evalued on the full set of low-order hybrid topologies (not inlcuding the pyramid). In this example, the temperature field is again, .. math:: :label: threed-L - + T = \frac{\lambda}{4} (cos(2 a \pi x) + cos(2 a \pi y) + cos(2 a \pi z)). Figure :numref:`laplace-hybrid` represents the MMS field for temperature @@ -843,7 +841,7 @@ verified by a simple algebraic calculation. .. code-block:: yaml actuator: - type: ActLineSimple + type: ActLineSimple search_method: stk_kdtree search_target_part: Unspecified-2-HEX @@ -898,7 +896,7 @@ conditions are given by :math:`\rho=1.0\textrm{ kg/m}^3` and :math:`U=2\textrm{ m/s}`. The aerodynamic properties of the fixed wing are given by the linear -lift coefficient +lift coefficient .. math:: :label: linear_cl @@ -968,13 +966,13 @@ work is an infinitesimally thin wing with a maximum chord (:math:`c_0`) of .. math:: :label: ew_sim_params - + \textrm{Span } b = \; & 10m \\ \textrm{Max chord } c_0 = \; & 1.0m \\ \textrm{Angle of attack } \alpha = \; & 7^{\circ} \\ U_{\infty} = \; & 10.0m/s \\ Reynolds \textrm{ number based on chord } = \; & 0.66M \\ - \textrm{Number of actuator points across span } = \; & 50 + \textrm{Number of actuator points across span } = \; & 50 The flow past the elliptic wing is simulated in a domain of size @@ -996,7 +994,7 @@ compared to simulations *a,b,c*. .. table:: - + +---------------+---------------------+----------------+------------------------+ | Case | :math:`\Delta x/c_0`|:math:`\Delta t`|:math:`\epsilon/\Delta` | +===============+=====================+================+========================+ @@ -1033,7 +1031,7 @@ drag per unit length along the blade. :numref:`ew_aoa` shows the comparison of the predicted angle of attack on the blade to the constant angle attack predicted by the lifting line theory. As expected, the agreement with the lifting line theory is much better near the -mid-span region compared to the wing tips. +mid-span region compared to the wing tips. .. _ew_cl: @@ -1088,24 +1086,24 @@ option, or 2) to use the standard open boundary condition in which the buoyancy reference value, rather than a single reference value. We test these open boundary condition options on a simplified stratified flow through a channel with slip walls. The -flow entering the domain is non-turbulent and uniformly 8 m/s. The temperature linearly varies from 300 K to 310 K from +flow entering the domain is non-turbulent and uniformly 8 m/s. The temperature linearly varies from 300 K to 310 K from the bottom to top of the channel with compatible, opposite-sign heat flux on the two walls to maintain this profile. -The Boussinesq buoyancy option is used, and the density is set constant to 1.17804 kg/m :math:`^3`. This density is -compatible with the reference pressure of 101325 Pa and a reference temperature of 300 K. The viscosity is set to +The Boussinesq buoyancy option is used, and the density is set constant to 1.17804 kg/m :math:`^3`. This density is +compatible with the reference pressure of 101325 Pa and a reference temperature of 300 K. The viscosity is set to 1.0e-5 Pa-s. *The flow should keep its inflow velocity and temperature profiles throughout the length of the domain*. -The domain is 3000 m long, 1000 m tall, and 20 m wide with 300 x 100 x 2 elements. The upper and lower boundaries +The domain is 3000 m long, 1000 m tall, and 20 m wide with 300 x 100 x 2 elements. The upper and lower boundaries are symmetry with the specified normal gradient of temperature option used such that the gradient matches the initial -temperature profile with its gradient of 0.01 K/m. Flow enters from the left and exits on the right. The remaining -boundaries are periodic. +temperature profile with its gradient of 0.01 K/m. Flow enters from the left and exits on the right. The remaining +boundaries are periodic. -We test the problem on three configurations: 1) using the standard open boundary condition, +We test the problem on three configurations: 1) using the standard open boundary condition, 2) using the global-mass-flow-rate-correction option, and 3) using the standard open boundary condition with a local moving-time-averaged reference temperature in the Boussinesq buoyancy term. Figure :numref:`stratified_outflow_ux1` shows the across-channel profile of outflow streamwise velocity. It is clear -that in configuration 1, the velocity is significantly distorted from the correct solution. Configurations 2 and 3 -remedy the problem. However, if we reduce the range of the x-axis, as shown in Figure :numref:`stratified_outflow_ux2`, +that in configuration 1, the velocity is significantly distorted from the correct solution. Configurations 2 and 3 +remedy the problem. However, if we reduce the range of the x-axis, as shown in Figure :numref:`stratified_outflow_ux2`, we see that configuration 3, the use of the standard open boundary condition with a local moving-time-averaged Boussinesq reference temperature, provides a superior solution in this case. In Figure, :numref:`stratified_outflow_T1`, we also see that configuration 1 significantly distorts the temperature from the correct solution. @@ -1139,12 +1137,12 @@ we also see that configuration 1 significantly distorts the temperature from the Outflow temperature profiles for the thermally stratified slip-channel flow. -We also verify that the global mass-flow-rate correction of configuration 2 is correcting the outflow mass flow rate +We also verify that the global mass-flow-rate correction of configuration 2 is correcting the outflow mass flow rate properly. The output from Nalu-Wind showing the correction is correct and is shown as follows: .. code-block:: c++ - + Mass Balance Review: Density accumulation: 0 Integrated inflow: -188486.0356751138 @@ -1158,18 +1156,18 @@ properly. The output from Nalu-Wind showing the correction is correct and is sh Specified Normal Temperature Gradient Boundary Condition -------------------------------------------------------- -The motivation for adding the ability to specify the boundary-normal temperature +The motivation for adding the ability to specify the boundary-normal temperature gradient is atmospheric boundary layer simulation in which the upper portion of the domain often contains a stably stratified layer with a temperature gradient that extends all the way to the upper boundary. The desire is for the simulation -to maintain that gradient throughout the simulation duration. +to maintain that gradient throughout the simulation duration. -Our test case is a laminar infinite channel with slip walls. In this case, the +Our test case is a laminar infinite channel with slip walls. In this case, the flow velocity is zero so the problem is simply a heat conduction through fluid. The density is fixed as constant, and there are no source terms including buoyancy. -This problem has an the analytical solution for the temperature profile across +This problem has an the analytical solution for the temperature profile across the channel: .. math:: @@ -1177,17 +1175,17 @@ the channel: T(t,z) = T(t_0,z_0) + \frac{-g_H-g_0}{H} \kappa_{eff} (t-t_0) + g_0 (z-z_0) + \frac{-g_H-g_0}{2H} (z-z_0)^2, -where :math:`t_0` is the initial time; :math:`z_0` is the height of the lower +where :math:`t_0` is the initial time; :math:`z_0` is the height of the lower channel wall; :math:`H` is the channel height; :math:`g_0` and :math:`g_H` are the wall-normal gradients of temperature at the lower and upper walls, respectively; :math:`\kappa_{eff}` is the effective thermal diffusivity; -and :math:`z` is the distance in the cross-channel direction. The sign of the +and :math:`z` is the distance in the cross-channel direction. The sign of the temperature gradients assumes that boundary normal points inward from the boundary. -For this solution to hold, the initial solution must be that of :eq:`T-slip-channel` +For this solution to hold, the initial solution must be that of :eq:`T-slip-channel` with :math:`t=t_0`. For all test cases, we use a domain that is 10 m x 10 m in the periodic (infinite) directions, -and 100 m in the cross-channel (z) direction. We specify a constant density of +and 100 m in the cross-channel (z) direction. We specify a constant density of 1 kg/m :math:`^3`, zero velocity, no buoyancy source term, a viscosity of 1 Pa-s, and a laminar Prandtl number of 1. No turbulence model is used. The value of :math:`T(t_0,z_0)` is 300 K. @@ -1198,7 +1196,7 @@ Simple Linear Temperature Profile: Equal and Opposite Specified Temperature Grad A simple verification test that is representative of a stable atmospheric capping inversion is to compute the simple thermal channel with equal and opposite specified -temperature gradients on each wall. By setting :math:`g_H = - g_0` in Equation +temperature gradients on each wall. By setting :math:`g_H = - g_0` in Equation :eq:`T-slip-channel`, we are left with .. math:: @@ -1212,7 +1210,7 @@ we set :math:`g_0 = 0.01` K/m and :math:`g_H = -0.01` K/m. We use a mesh that 2 elements wide in the periodic directions and 20 elements across the channel. We simulate a long time period of 25,000 s. Figure :numref:`T_gradBC_linear` -shows that the computed and analytical solutions agree. +shows that the computed and analytical solutions agree. .. _T_gradBC_linear: @@ -1220,14 +1218,14 @@ shows that the computed and analytical solutions agree. :width: 500px :align: center - The analytical (black solid) and computed (red dashed) temperature profile from + The analytical (black solid) and computed (red dashed) temperature profile from the case with :math:`g_H = -g_0` at :math:`t =` 25,000 s. Parabolic Temperature Profile: Equal Specified Temperature Gradients ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -Next, we verify the specified normal temperature gradient boundary condition +Next, we verify the specified normal temperature gradient boundary condition option by computing the simple thermal channel with equal specified temperature gradients, which yields the full time-dependent solution of Equation :eq:`T-slip-channel`. Here, we set :math:`g_0 = g_H = 0.01` K/m. @@ -1244,5 +1242,5 @@ coarser meshes. :width: 500px :align: center - The analytical (black solid) and computed (colored) temperature profile from + The analytical (black solid) and computed (colored) temperature profile from the case with :math:`g_H = g_0` at :math:`t =` 25,000 s. diff --git a/include/FieldRegistry.h b/include/FieldRegistry.h index 5067f506d..0eef91024 100644 --- a/include/FieldRegistry.h +++ b/include/FieldRegistry.h @@ -41,7 +41,7 @@ class FieldRegistry break; } default: - throw std::runtime_error("Unsupported number of reference states"); + throw std::runtime_error("Unsupported number of reference states for field " + name + " with " + std::to_string(numStates) + " states and " + std::to_string(numDim) + " dims."); } break; } @@ -56,7 +56,7 @@ class FieldRegistry break; } default: - throw std::runtime_error("Unsupported number of reference states"); + throw std::runtime_error("Unsupported number of reference states for field " + name + " with " + std::to_string(numStates) + " states and " + std::to_string(numDim) + " dims."); } break; } diff --git a/unit_tests/algorithms/CMakeLists.txt b/unit_tests/algorithms/CMakeLists.txt index 0410b8884..b2927cfb9 100644 --- a/unit_tests/algorithms/CMakeLists.txt +++ b/unit_tests/algorithms/CMakeLists.txt @@ -2,4 +2,5 @@ target_sources(${utest_ex_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestAlgorithm.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestLESAlgorithms.C ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestSSTAlgorithms.C + ${CMAKE_CURRENT_SOURCE_DIR}/UnitTestMomentumBoussinesqSrcNodeSuppAlg.C ) diff --git a/unit_tests/algorithms/UnitTestMomentumBoussinesqSrcNodeSuppAlg.C b/unit_tests/algorithms/UnitTestMomentumBoussinesqSrcNodeSuppAlg.C index b0dcccfa6..535a3e5c2 100644 --- a/unit_tests/algorithms/UnitTestMomentumBoussinesqSrcNodeSuppAlg.C +++ b/unit_tests/algorithms/UnitTestMomentumBoussinesqSrcNodeSuppAlg.C @@ -17,7 +17,6 @@ #include "UnitTestRealm.h" #include "Realm.h" #include "SolutionOptions.h" -#include "MomentumBoussinesqSrcNodeSuppAlg.h" #include "MomentumBoussinesqRASrcNodeSuppAlg.h" #include "MovingAveragePostProcessor.h" #include "TimeIntegrator.h" @@ -26,58 +25,7 @@ #include #include -TEST(MomentumBoussinesqSrcNodeSuppAlg, single_value) -{ - NodeSuppHelper helper; - auto& meta = helper.realm.meta_data(); - - auto& dnv = - meta.declare_field(stk::topology::NODE_RANK, "dual_nodal_volume"); - stk::mesh::put_field_on_mesh(dnv, meta.universal_part(), nullptr); - - auto& temperature = - meta.declare_field(stk::topology::NODE_RANK, "temperature"); - stk::mesh::put_field_on_mesh(temperature, meta.universal_part(), nullptr); - - meta.commit(); - - stk::mesh::Entity node = helper.make_one_node_mesh(); - - const double dnv_value = 0.125; - *stk::mesh::field_data(dnv, node) = dnv_value; - - const double temperature_value = 305; - *stk::mesh::field_data(temperature, node) = temperature_value; - - auto& solnOpts = *helper.realm.solutionOptions_; - - const double tRef = 300; - solnOpts.referenceTemperature_ = tRef; - - const double rhoRef = 1.0; - solnOpts.referenceDensity_ = rhoRef; - - const double beta = 1.0 / tRef; - solnOpts.thermalExpansionCoeff_ = beta; - - std::vector gravity = {-5, 6, 7}; - solnOpts.gravity_ = gravity; - - double coeff = -beta * rhoRef * (temperature_value - tRef) * dnv_value; - double expected_rhs[3] = { - coeff * gravity[0], coeff * gravity[1], coeff * gravity[2]}; - - double rhs[3] = {0, 0, 0}; - auto boussinesqAlg = - sierra::nalu::MomentumBoussinesqSrcNodeSuppAlg(helper.realm); - boussinesqAlg.node_execute(nullptr, rhs, node); - - for (int d = 0; d < 3; ++d) { - EXPECT_DOUBLE_EQ(rhs[d], expected_rhs[d]); - } -} - -TEST(MomentumBoussinesqRASrcNodeSuppAlg, single_value) +TEST(TestMomentumBoussinesqRASrcNodeSuppAlg, single_value) { NodeSuppHelper helper; auto& meta = helper.realm.meta_data(); diff --git a/unit_tests/algorithms/UnitTestRodiAlgorithm.C b/unit_tests/algorithms/UnitTestRodiAlgorithm.C deleted file mode 100644 index efb9b805e..000000000 --- a/unit_tests/algorithms/UnitTestRodiAlgorithm.C +++ /dev/null @@ -1,52 +0,0 @@ -// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC -// (NTESS), National Renewable Energy Laboratory, University of Texas Austin, -// Northwest Research Associates. Under the terms of Contract DE-NA0003525 -// with NTESS, the U.S. Government retains certain rights in this software. -// -// This software is released under the BSD 3-clause license. See LICENSE file -// for more details. -// - -#include "UnitTestAlgorithm.h" -#include "UnitTestKokkosUtils.h" -#include "UnitTestFieldUtils.h" -#include "UnitTestAlgorithmUtils.h" - -#include "SolutionOptions.h" -#include "TurbKineticEnergyRodiNodeSourceSuppAlg.h" - -TEST_F(TestTurbulenceAlgorithm, turbkineticenergyrodinodesourcesuppalg) -{ - sierra::nalu::Realm& realm = this->create_realm(); - - const int nprocs = this->bulk().parallel_size(); - std::string meshSpec = "generated:1x1x" + std::to_string(nprocs); - fill_mesh_and_init_fields(meshSpec); - - // set solution options - realm.solutionOptions_->gravity_.resize(3); - realm.solutionOptions_->gravity_[0] = 10.0; - realm.solutionOptions_->gravity_[1] = -10.0; - realm.solutionOptions_->gravity_[2] = 5.0; - realm.solutionOptions_->turbPrMap_["enthalpy"] = 0.60; - realm.solutionOptions_->thermalExpansionCoeff_ = 3.0e-3; - - // Nodal execute - auto& bulk = this->bulk(); - unit_test_algorithm_utils::TestSupplementalAlgorithmDriver assembleSuppAlgs( - bulk); - std::unique_ptr suppalg( - new sierra::nalu::TurbKineticEnergyRodiNodeSourceSuppAlg(realm)); - assembleSuppAlgs.activeSuppAlgs_.push_back(suppalg.get()); - assembleSuppAlgs.nodal_execute(); - - // Perform tests - const double tol = 1e-14; - const double lhs_norm = assembleSuppAlgs.get_lhs_norm(); - const double rhs_norm = assembleSuppAlgs.get_rhs_norm(); - const double lhs_gold_norm = 0.0; - const double rhs_gold_norm = - nprocs == 1 ? 0.00030133929246584567 : 0.0002760523778561557; - EXPECT_NEAR(lhs_norm, lhs_gold_norm, tol); - EXPECT_NEAR(rhs_norm, rhs_gold_norm, tol); -}