Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/Problem/InterfaceAdvection.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Problem configuration for the compressible flow equations
9 : : \details This file defines a Problem policy class for the compressible flow
10 : : equations, defined in PDE/MultiMat/MultiMat.h. See PDE/MultiMat/Problem.h
11 : : for general requirements on Problem policy classes for MultiMat.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include "InterfaceAdvection.hpp"
16 : : #include "Inciter/InputDeck/InputDeck.hpp"
17 : : #include "EoS/EoS.hpp"
18 : : #include "MultiMat/MultiMatIndexing.hpp"
19 : :
20 : : using inciter::MultiMatProblemInterfaceAdvection;
21 : :
22 : : tk::InitializeFn::result_type
23 : 419586 : MultiMatProblemInterfaceAdvection::initialize( ncomp_t system,
24 : : ncomp_t ncomp,
25 : : tk::real x,
26 : : tk::real y,
27 : : tk::real /*z*/,
28 : : tk::real t )
29 : : // *****************************************************************************
30 : : //! Evaluate analytical solution at (x,y,z,t) for all components
31 : : //! \param[in] system Equation system index, i.e., which compressible
32 : : //! flow equation system we operate on among the systems of PDEs
33 : : //! \param[in] ncomp Number of scalar components in this PDE system
34 : : //! \param[in] x X coordinate where to evaluate the solution
35 : : //! \param[in] y Y coordinate where to evaluate the solution
36 : : //! \param[in] t Time where to evaluate the solution
37 : : //! \return Values of all components evaluated at (x)
38 : : //! \note The function signature must follow tk::InitializeFn
39 : : // *****************************************************************************
40 : : {
41 : : auto nmat =
42 : 419586 : g_inputdeck.get< tag::param, eq, tag::nmat >()[system];
43 : :
44 : : // see also Control/Inciter/InputDeck/Grammar.hpp
45 : : Assert( ncomp == 3*nmat+3, "Incorrect number of components in multi-material "
46 : : "system" );
47 : :
48 : 419586 : std::vector< tk::real > s( ncomp, 0.0 );
49 : : auto u = std::sqrt(50.0);
50 : : auto v = std::sqrt(50.0);
51 : : auto w = 0.0;
52 : : auto alphamin = 1.0e-12;
53 : :
54 : : // center of the cylinder
55 : 419586 : auto x0 = 0.45 + u*t;
56 : : auto y0 = 0.45 + v*t;
57 : :
58 : : // radii of the material-rings
59 [ + - ][ - - ]: 419586 : std::vector< tk::real > r0(nmat, 0.0);
60 : 419586 : r0[nmat-1] = 0.0;
61 : 419586 : r0[nmat-2] = 0.1;
62 : 419586 : r0[0] = 0.35;
63 [ - + ]: 419586 : for (std::size_t k=1; k<nmat-2; ++k)
64 [ - - ]: 0 : r0[k] = r0[k-1] - (r0[0]-r0[nmat-2])
65 [ - - ]: 0 : /(std::max( 1.0, static_cast<tk::real>(nmat-2)) );
66 : :
67 [ + + ]: 1678344 : for (std::size_t k=0; k<nmat; ++k)
68 : 1258758 : s[volfracIdx(nmat, k)] = alphamin;
69 : :
70 : : // interface location
71 : 419586 : auto r = std::sqrt( (x-x0)*(x-x0) + (y-y0)*(y-y0) );
72 : : bool is_mat(false);
73 [ + + ]: 1258758 : for (std::size_t k=0; k<nmat-1; ++k)
74 : : {
75 [ + + ][ + + ]: 839172 : if (r<r0[k] && r>=r0[k+1])
76 : : {
77 : 132702 : s[volfracIdx(nmat, k)] = 1.0 - static_cast<tk::real>(nmat-1)*alphamin;
78 : : is_mat = true;
79 : : }
80 : : }
81 [ + + ]: 419586 : if (!is_mat)
82 : : {
83 : 286884 : s[volfracIdx(nmat, nmat-1)] = 1.0 - static_cast<tk::real>(nmat-1)*alphamin;
84 : : }
85 : :
86 : : auto rhob = 0.0;
87 [ + + ]: 1678344 : for (std::size_t k=0; k<nmat; ++k)
88 : : {
89 : 1258758 : auto rhok = eos_density< eq >( system, 1.0e5, 300.0, k );
90 : 1258758 : s[densityIdx(nmat, k)] = s[volfracIdx(nmat, k)] * rhok;
91 : 1258758 : s[energyIdx(nmat, k)] = s[volfracIdx(nmat, k)]
92 : 1258758 : * eos_totalenergy< eq >( system, rhok, u, v, w, 1.0e5, k );
93 : 1258758 : rhob += s[densityIdx(nmat, k)];
94 : : }
95 [ + - ]: 419586 : s[momentumIdx(nmat, 0)] = rhob * u;
96 : 419586 : s[momentumIdx(nmat, 1)] = rhob * v;
97 [ + - ]: 419586 : s[momentumIdx(nmat, 2)] = rhob * w;
98 : :
99 : 419586 : return s;
100 : : }
101 : :
102 : : std::vector< std::string >
103 : 2 : MultiMatProblemInterfaceAdvection::names( ncomp_t )
104 : : // *****************************************************************************
105 : : // Return names of integral variables to be output to diagnostics file
106 : : //! \return Vector of strings labelling integral variables output
107 : : // *****************************************************************************
108 : : {
109 [ + - ][ + - ]: 12 : return { "r", "ru", "rv", "rw", "re" };
[ + - ][ + - ]
[ + - ][ - + ]
[ + + ][ - + ]
[ - - ][ - - ]
110 : : }
|
