Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Transport/Problem/SlotCyl.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 transport equations
9 : : \details This file defines a Problem policy class for the transport
10 : : equations, defined in PDE/Transport/CGTransport.h implementing
11 : : node-centered continuous Galerkin (CG) and PDE/Transport/DGTransport.h
12 : : implementing cell-centered discontinuous Galerkin (DG) discretizations.
13 : : See PDE/Transport/Problem.h for general requirements on Problem policy
14 : : classes for cg::Transport and dg::Transport.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include "SlotCyl.hpp"
19 : : #include "SystemComponents.hpp"
20 : : #include "Inciter/InputDeck/InputDeck.hpp"
21 : :
22 : : namespace inciter {
23 : :
24 : : extern ctr::InputDeck g_inputdeck;
25 : :
26 : : } // ::inciter
27 : :
28 : : using inciter::TransportProblemSlotCyl;
29 : :
30 : : std::vector< tk::real >
31 : 15037973 : TransportProblemSlotCyl::initialize( ncomp_t, ncomp_t ncomp, tk::real x,
32 : : tk::real y, tk::real, tk::real t )
33 : : // *****************************************************************************
34 : : // Evaluate analytical solution at (x,y,t) for all components
35 : : //! \param[in] ncomp Number of components in this transport equation system
36 : : //! \param[in] x X coordinate where to evaluate the solution
37 : : //! \param[in] y Y coordinate where to evaluate the solution
38 : : //! \param[in] t Time where to evaluate the solution
39 : : //! \return Values of all components evaluated at (x,y,t)
40 : : // *****************************************************************************
41 : : {
42 : : using std::sin; using std::cos;
43 : :
44 : 15037973 : std::vector< tk::real > s( ncomp, 0.0 );
45 [ + + ]: 30075946 : for (ncomp_t c=0; c<ncomp; ++c) {
46 : 15037973 : auto T = t + 2.0*M_PI/static_cast<tk::real>(ncomp)*static_cast<tk::real>(c);
47 : : const tk::real R0 = 0.15;
48 : :
49 : : // center of the cone
50 : : tk::real x0 = 0.5;
51 : : tk::real y0 = 0.25;
52 : : tk::real r = std::sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
53 : 15037973 : tk::real kx = 0.5 + r*sin( T );
54 : 15037973 : tk::real ky = 0.5 - r*cos( T );
55 : :
56 : : // center of the hump
57 : : x0 = 0.25;
58 : : y0 = 0.5;
59 : : r = std::sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
60 : 15037973 : tk::real hx = 0.5 + r*sin( T-M_PI/2.0 ),
61 : 15037973 : hy = 0.5 - r*cos( T-M_PI/2.0 );
62 : :
63 : : // center of the slotted cylinder
64 : : x0 = 0.5;
65 : : y0 = 0.75;
66 : : r = std::sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
67 : 15037973 : tk::real cx = 0.5 + r*sin( T+M_PI ),
68 : 15037973 : cy = 0.5 - r*cos( T+M_PI );
69 : :
70 : : // end points of the cylinder slot
71 : 15037973 : tk::real i1x = 0.525, i1y = cy - r*cos( std::asin(0.025/r) ),
72 : : i2x = 0.525, i2y = 0.8,
73 : : i3x = 0.475, i3y = 0.8;
74 : :
75 : : // rotate end points of cylinder slot
76 : 15037973 : tk::real ri1x = 0.5 + cos(T)*(i1x-0.5) - sin(T)*(i1y-0.5),
77 : 15037973 : ri1y = 0.5 + sin(T)*(i1x-0.5) + cos(T)*(i1y-0.5),
78 : 15037973 : ri2x = 0.5 + cos(T)*(i2x-0.5) - sin(T)*(i2y-0.5),
79 : 15037973 : ri2y = 0.5 + sin(T)*(i2x-0.5) + cos(T)*(i2y-0.5),
80 : 15037973 : ri3x = 0.5 + cos(T)*(i3x-0.5) - sin(T)*(i3y-0.5),
81 : 15037973 : ri3y = 0.5 + sin(T)*(i3x-0.5) + cos(T)*(i3y-0.5);
82 : :
83 : : // direction of slot sides
84 : 15037973 : tk::real v1x = ri2x-ri1x, v1y = ri2y-ri1y,
85 : 15037973 : v2x = ri3x-ri2x, v2y = ri3y-ri2y;
86 : :
87 : : // lengths of direction of slot sides vectors
88 : 15037973 : tk::real v1 = std::sqrt(v1x*v1x + v1y*v1y),
89 : 15037973 : v2 = std::sqrt(v2x*v2x + v2y*v2y);
90 : :
91 : : // cone
92 : 15037973 : r = std::sqrt((x-kx)*(x-kx) + (y-ky)*(y-ky)) / R0;
93 [ + + ]: 15037973 : if (r<1.0) s[c] = 0.6*(1.0-r);
94 : :
95 : : // hump
96 : 15037973 : r = std::sqrt((x-hx)*(x-hx) + (y-hy)*(y-hy)) / R0;
97 [ + + ][ - + ]: 15037973 : if (r<1.0) s[c] = 0.2*(1.0+cos(M_PI*std::min(r,1.0)));
98 : :
99 : : // cylinder
100 : 15037973 : r = std::sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy)) / R0;
101 : : const std::array< tk::real, 2 > r1{{ v1x, v1y }},
102 : 15037973 : r2{{ x-ri1x, y-ri1y }};
103 : 15037973 : const auto d1 = (r1[0]*r2[1] - r2[0]*r1[1]) / v1;
104 : : const std::array< tk::real, 2 > r3{{ v2x, v2y }},
105 : 15037973 : r4{{ x-ri2x, y-ri2y }};
106 : 15037973 : const auto d2 = (r3[0]*r4[1] - r4[0]*r3[1]) / v2;
107 [ + + ][ + + ]: 15037973 : if (r<1.0 && (d1>0.05 || d1<0.0 || d2<0.0)) s[c] = 0.6;
[ + + ][ + + ]
108 : : }
109 : :
110 : 15037973 : return s;
111 : : }
112 : :
113 : : std::vector< std::array< tk::real, 3 > >
114 : 128197793 : TransportProblemSlotCyl::prescribedVelocity( ncomp_t, ncomp_t ncomp,
115 : : tk::real x, tk::real y, tk::real,
116 : : tk::real )
117 : : // *****************************************************************************
118 : : // Assign prescribed shear velocity at a point
119 : : //! \param[in] ncomp Number of components in this transport equation
120 : : //! \param[in] x X coordinate at which to assign velocity
121 : : //! \param[in] y y coordinate at which to assign velocity
122 : : //! \return Velocity assigned to all vertices of a tetrehedron, size:
123 : : //! ncomp * ndim = [ncomp][3]
124 : : // *****************************************************************************
125 : : {
126 : 128197793 : std::vector< std::array< tk::real, 3 > > vel( ncomp );
127 : :
128 [ + + ]: 256395586 : for (ncomp_t c=0; c<ncomp; ++c)
129 : 128197793 : vel[c] = {{ 0.5-y, x-0.5, 0.0 }};
130 : :
131 : 128197793 : return vel;
132 : : }
|
