Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Zalesak.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Zalesak, FCT limiting for edge-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Vector.hpp"
14 : : #include "Around.hpp"
15 : : #include "DerivedData.hpp"
16 : : #include "EOS.hpp"
17 : : #include "Zalesak.hpp"
18 : : #include "Problems.hpp"
19 : : #include "InciterConfig.hpp"
20 : :
21 : : namespace inciter {
22 : :
23 : : extern ctr::Config g_cfg;
24 : :
25 : : } // ::inciter
26 : :
27 : : namespace zalesak {
28 : :
29 : : using inciter::g_cfg;
30 : :
31 : : static void
32 : 9794900 : advedge( const tk::real supint[],
33 : : const tk::Fields& U,
34 : : const std::array< std::vector< tk::real >, 3 >& coord,
35 : : tk::real t,
36 : : tk::real dt,
37 : : const std::vector< tk::real >& tp,
38 : : const std::vector< tk::real >& dtp,
39 : : std::size_t p,
40 : : std::size_t q,
41 : : tk::real f[],
42 : : const std::function< std::vector< tk::real >
43 : : ( tk::real, tk::real, tk::real, tk::real ) >& src )
44 : : // *****************************************************************************
45 : : //! Compute advection fluxes on a single edge
46 : : //! \param[in] supint Edge integral
47 : : //! \param[in] U Solution vector to read conserved variables from
48 : : //! \param[in] coord Mesh node coordinates
49 : : //! \param[in] t Physical time
50 : : //! \param[in] dt Physical time step size
51 : : //! \param[in] tp Phisical time step size for each mesh node (if steady state)
52 : : //! \param[in] dtp Time step size for each mesh node (if steady state)
53 : : //! \param[in] p Left node index of edge
54 : : //! \param[in] q Right node index of edge
55 : : //! \param[in,out] f Flux computed
56 : : //! \param[in] src Function to call to evaluate a problem-sepcific source term
57 : : // *****************************************************************************
58 : : {
59 : 9794900 : const auto steady = g_cfg.get< tag::steady >();
60 : : const auto ncomp = U.nprop();
61 : : const auto& x = coord[0];
62 : : const auto& y = coord[1];
63 : : const auto& z = coord[2];
64 : :
65 : : // edge vector
66 : 9794900 : auto dx = x[p] - x[q];
67 : 9794900 : auto dy = y[p] - y[q];
68 : 9794900 : auto dz = z[p] - z[q];
69 : 9794900 : auto dl = dx*dx + dy*dy + dz*dz;
70 : 9794900 : dx /= dl;
71 : 9794900 : dy /= dl;
72 [ + + ]: 9794900 : dz /= dl;
73 : :
74 : : // left state
75 : 9794900 : auto rL = U(p,0);
76 : 9794900 : auto ruL = U(p,1);
77 : 9794900 : auto rvL = U(p,2);
78 : 9794900 : auto rwL = U(p,3);
79 : 9794900 : auto reL = U(p,4);
80 [ + + ]: 9794900 : auto pL = eos::pressure( reL - 0.5*(ruL*ruL + rvL*rvL + rwL*rwL)/rL );
81 : 9794900 : auto dnL = (ruL*dx + rvL*dy + rwL*dz)/rL;
82 : :
83 : : // right state
84 : 9794900 : auto rR = U(q,0);
85 : 9794900 : auto ruR = U(q,1);
86 : 9794900 : auto rvR = U(q,2);
87 : 9794900 : auto rwR = U(q,3);
88 : 9794900 : auto reR = U(q,4);
89 : 9794900 : auto pR = eos::pressure( reR - 0.5*(ruR*ruR + rvR*rvR + rwR*rwR)/rR );
90 : 9794900 : auto dnR = (ruR*dx + rvR*dy + rwR*dz)/rR;
91 : :
92 : 9794900 : auto nx = supint[0];
93 : 9794900 : auto ny = supint[1];
94 : 9794900 : auto nz = supint[2];
95 : :
96 : : #if defined(__clang__)
97 : : #pragma clang diagnostic push
98 : : #pragma clang diagnostic ignored "-Wvla"
99 : : #pragma clang diagnostic ignored "-Wvla-extension"
100 : : #elif defined(STRICT_GNUC)
101 : : #pragma GCC diagnostic push
102 : : #pragma GCC diagnostic ignored "-Wvla"
103 : : #endif
104 : :
105 : : // Taylor-Galerkin first half step
106 : :
107 [ + + ]: 9794900 : if (steady) dt = (dtp[p] + dtp[q])/2.0;
108 : :
109 : 9794900 : tk::real ue[ncomp];
110 : :
111 : : // flow
112 : 9794900 : auto dp = pL - pR;
113 : 9794900 : ue[0] = 0.5*(rL + rR - dt*(rL*dnL - rR*dnR));
114 : 9794900 : ue[1] = 0.5*(ruL + ruR - dt*(ruL*dnL - ruR*dnR + dp*dx));
115 : 9794900 : ue[2] = 0.5*(rvL + rvR - dt*(rvL*dnL - rvR*dnR + dp*dy));
116 : 9794900 : ue[3] = 0.5*(rwL + rwR - dt*(rwL*dnL - rwR*dnR + dp*dz));
117 : 9794900 : ue[4] = 0.5*(reL + reR - dt*((reL+pL)*dnL - (reR+pR)*dnR));
118 : : // scalar
119 [ + + ]: 10192880 : for (std::size_t c=5; c<ncomp; ++c) {
120 : 397980 : ue[c] = 0.5*(U(p,c) + U(q,c) - dt*(U(p,c)*dnL - U(q,c)*dnR));
121 : : }
122 : :
123 : : // source
124 [ + + ]: 9794900 : if (src) {
125 [ - + ]: 423100 : if (steady) t = (tp[p] + tp[q])/2.0;
126 : 423100 : auto coef = dt/4.0;
127 [ - + ]: 423100 : auto sL = src( x[p], y[p], z[p], t );
128 [ - + ]: 423100 : auto sR = src( x[q], y[q], z[q], t );
129 : : // flow + scalar
130 [ + + ]: 2936580 : for (std::size_t c=0; c<ncomp; ++c) {
131 : 2513480 : ue[c] += coef*(sL[c] + sR[c]);
132 : : }
133 : : }
134 : :
135 : : // Taylor-Galerkin second half step
136 : :
137 : 9794900 : auto rh = ue[0];
138 : 9794900 : auto ruh = ue[1];
139 : 9794900 : auto rvh = ue[2];
140 : 9794900 : auto rwh = ue[3];
141 : 9794900 : auto reh = ue[4];
142 : 9794900 : auto ph = eos::pressure( reh - 0.5*(ruh*ruh + rvh*rvh + rwh*rwh)/rh );
143 : 9794900 : auto vn = (ruh*nx + rvh*ny + rwh*nz)/rh;
144 : :
145 : : // flow
146 : 9794900 : f[0] = 2.0*rh*vn;
147 : 9794900 : f[1] = 2.0*(ruh*vn + ph*nx);
148 : 9794900 : f[2] = 2.0*(rvh*vn + ph*ny);
149 : 9794900 : f[3] = 2.0*(rwh*vn + ph*nz);
150 : 9794900 : f[4] = 2.0*(reh + ph)*vn;
151 : : // scalar
152 [ + + ]: 10192880 : for (std::size_t c=5; c<ncomp; ++c) {
153 : 397980 : f[c] = 2.0*ue[c]*vn;
154 : : }
155 : :
156 : : // source
157 [ + + ]: 9794900 : if (src) {
158 : 423100 : auto coef = -5.0/3.0*supint[3];
159 : 423100 : auto xe = (x[p] + x[q])/2.0;
160 : 423100 : auto ye = (y[p] + y[q])/2.0;
161 : 423100 : auto ze = (z[p] + z[q])/2.0;
162 : 423100 : auto se = src( xe, ye, ze, t+dt/2.0 );
163 : : // flow + scalar
164 [ + + ]: 2936580 : for (std::size_t c=0; c<ncomp; ++c) {
165 : 2513480 : f[ncomp+c] = coef*se[c];
166 : : }
167 : : }
168 : :
169 : : // artificial viscosity
170 : :
171 : 9794900 : const auto stab2 = g_cfg.get< tag::stab2 >();
172 [ + + ]: 9794900 : if (!stab2) return;
173 : :
174 : 5838840 : auto stab2coef = g_cfg.get< tag::stab2coef >();
175 : 5838840 : auto vnL = (ruL*nx + rvL*ny + rwL*nz)/rL;
176 [ + + ]: 5838840 : auto vnR = (ruR*nx + rvR*ny + rwR*nz)/rR;
177 : : auto len = tk::length( nx, ny, nz );
178 [ + + ][ - + ]: 5839250 : auto cL = eos::soundspeed( std::max(rL,1.0e-8), std::max(pL,0.0) );
[ + + ]
179 [ + + ][ - + ]: 5839413 : auto cR = eos::soundspeed( std::max(rR,1.0e-8), std::max(pR,0.0) );
[ + + ]
180 : 5838840 : auto sl = std::abs(vnL) + cL*len;
181 [ + + ]: 5838840 : auto sr = std::abs(vnR) + cR*len;
182 : 5838840 : auto fw = stab2coef * std::max( sl, sr );
183 : :
184 : : // flow
185 : 5838840 : f[0] -= fw*(rL - rR);
186 : 5838840 : f[1] -= fw*(ruL - ruR);
187 : 5838840 : f[2] -= fw*(rvL - rvR);
188 : 5838840 : f[3] -= fw*(rwL - rwR);
189 : 5838840 : f[4] -= fw*(reL - reR);
190 : : // scalar
191 [ - + ]: 5838840 : for (std::size_t c=5; c<ncomp; ++c) {
192 : 0 : f[c] -= fw*(U(p,c) - U(q,c));
193 : : }
194 : :
195 : : #if defined(__clang__)
196 : : #pragma clang diagnostic pop
197 : : #elif defined(STRICT_GNUC)
198 : : #pragma GCC diagnostic pop
199 : : #endif
200 [ + + ]: 9794900 : }
201 : :
202 : : static void
203 : 5365 : advdom( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
204 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
205 : : const std::array< std::vector< tk::real >, 3 >& coord,
206 : : tk::real t,
207 : : tk::real dt,
208 : : const std::vector< tk::real >& tp,
209 : : const std::vector< tk::real >& dtp,
210 : : const tk::Fields& U,
211 : : // cppcheck-suppress constParameter
212 : : tk::Fields& R )
213 : : // *****************************************************************************
214 : : //! Compute domain-edge integrals for advection
215 : : //! \param[in] dsupedge Domain superedges
216 : : //! \param[in] dsupint Domain superedge integrals
217 : : //! \param[in] coord Mesh node coordinates
218 : : //! \param[in] t Physical time
219 : : //! \param[in] dt Physical time step size
220 : : //! \param[in] tp Phisical time step size for each mesh node (if steady state)
221 : : //! \param[in] dtp Time step size for each mesh node (if steady state)
222 : : //! \param[in] U Solution vector at recent time step
223 : : //! \param[in,out] R Right-hand side vector
224 : : // *****************************************************************************
225 : : {
226 : : auto ncomp = U.nprop();
227 : 5365 : auto src = problems::SRC();
228 : :
229 : : #if defined(__clang__)
230 : : #pragma clang diagnostic push
231 : : #pragma clang diagnostic ignored "-Wvla"
232 : : #pragma clang diagnostic ignored "-Wvla-extension"
233 : : #elif defined(STRICT_GNUC)
234 : : #pragma GCC diagnostic push
235 : : #pragma GCC diagnostic ignored "-Wvla"
236 : : #endif
237 : :
238 : : // domain edge contributions: tetrahedron superedges
239 [ + + ]: 988400 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
240 : 983035 : const auto N = dsupedge[0].data() + e*4;
241 : 983035 : tk::real f[6][ncomp*2];
242 : : const auto d = dsupint[0].data();
243 [ + - ]: 983035 : advedge( d+(e*6+0)*4, U, coord, t, dt, tp, dtp, N[0], N[1], f[0], src );
244 [ + - ]: 983035 : advedge( d+(e*6+1)*4, U, coord, t, dt, tp, dtp, N[1], N[2], f[1], src );
245 [ + - ]: 983035 : advedge( d+(e*6+2)*4, U, coord, t, dt, tp, dtp, N[2], N[0], f[2], src );
246 [ + - ]: 983035 : advedge( d+(e*6+3)*4, U, coord, t, dt, tp, dtp, N[0], N[3], f[3], src );
247 [ + - ]: 983035 : advedge( d+(e*6+4)*4, U, coord, t, dt, tp, dtp, N[1], N[3], f[4], src );
248 [ + - ]: 983035 : advedge( d+(e*6+5)*4, U, coord, t, dt, tp, dtp, N[2], N[3], f[5], src );
249 [ + + ]: 5930990 : for (std::size_t c=0; c<ncomp; ++c) {
250 [ + + ]: 4947955 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
251 : 4947955 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
252 : 4947955 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
253 [ + + ]: 4947955 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
254 [ + + ]: 4947955 : if (src) {
255 : 206980 : auto nc = ncomp + c;
256 : 206980 : R(N[0],c) += f[0][nc] + f[2][nc] + f[3][nc];
257 : 206980 : R(N[1],c) += f[0][nc] + f[1][nc] + f[4][nc];
258 : 206980 : R(N[2],c) += f[1][nc] + f[2][nc] + f[5][nc];
259 : 206980 : R(N[3],c) += f[3][nc] + f[4][nc] + f[5][nc];
260 : : }
261 : : }
262 [ + - ]: 983035 : }
263 : :
264 : : // domain edge contributions: triangle superedges
265 [ + + ]: 545395 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
266 : 540030 : const auto N = dsupedge[1].data() + e*3;
267 : 540030 : tk::real f[3][ncomp*2];
268 : : const auto d = dsupint[1].data();
269 [ + - ]: 540030 : advedge( d+(e*3+0)*4, U, coord, t, dt, tp, dtp, N[0], N[1], f[0], src );
270 [ + - ]: 540030 : advedge( d+(e*3+1)*4, U, coord, t, dt, tp, dtp, N[1], N[2], f[1], src );
271 [ + - ]: 540030 : advedge( d+(e*3+2)*4, U, coord, t, dt, tp, dtp, N[2], N[0], f[2], src );
272 [ + + ]: 3271440 : for (std::size_t c=0; c<ncomp; ++c) {
273 [ + + ]: 2731410 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
274 : 2731410 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
275 [ + + ]: 2731410 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
276 [ + + ]: 2731410 : if (src) {
277 : 198110 : auto nc = ncomp + c;
278 : 198110 : R(N[0],c) += f[0][nc] + f[2][nc];
279 : 198110 : R(N[1],c) += f[0][nc] + f[1][nc];
280 : 198110 : R(N[2],c) += f[1][nc] + f[2][nc];
281 : : }
282 : : }
283 [ + - ]: 540030 : }
284 : :
285 : : // domain edge contributions: edges
286 [ + + ]: 2281965 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
287 : 2276600 : const auto N = dsupedge[2].data() + e*2;
288 : 2276600 : tk::real f[ncomp*2];
289 : : const auto d = dsupint[2].data();
290 [ + - ]: 2276600 : advedge( d+e*4, U, coord, t, dt, tp, dtp, N[0], N[1], f, src );
291 [ + + ]: 13767120 : for (std::size_t c=0; c<ncomp; ++c) {
292 [ + + ]: 11490520 : R(N[0],c) -= f[c];
293 [ + + ]: 11490520 : R(N[1],c) += f[c];
294 [ + + ]: 11490520 : if (src) {
295 : 677270 : auto nc = ncomp + c;
296 : 677270 : R(N[0],c) += f[nc];
297 : 677270 : R(N[1],c) += f[nc];
298 : : }
299 : : }
300 [ + - ]: 2276600 : }
301 : :
302 : : #if defined(__clang__)
303 : : #pragma clang diagnostic pop
304 : : #elif defined(STRICT_GNUC)
305 : : #pragma GCC diagnostic pop
306 : : #endif
307 : 5365 : }
308 : :
309 : : static void
310 : 5365 : advbnd( const std::vector< std::size_t >& triinpoel,
311 : : const std::array< std::vector< tk::real >, 3 >& coord,
312 : : const std::vector< std::uint8_t >& besym,
313 : : const tk::Fields& U,
314 : : tk::Fields& R )
315 : : // *****************************************************************************
316 : : //! Compute boundary integrals for advection
317 : : //! \param[in] triinpoel Boundary face connectivity
318 : : //! \param[in] coord Mesh node coordinates
319 : : //! \param[in] besym Boundary element symmetry BC flags
320 : : //! \param[in] U Solution vector at recent time step
321 : : //! \param[in,out] R Right-hand side vector
322 : : // *****************************************************************************
323 : : {
324 : : auto ncomp = U.nprop();
325 : :
326 : : const auto& x = coord[0];
327 : : const auto& y = coord[1];
328 : : const auto& z = coord[2];
329 : :
330 : : #if defined(__clang__)
331 : : #pragma clang diagnostic push
332 : : #pragma clang diagnostic ignored "-Wvla"
333 : : #pragma clang diagnostic ignored "-Wvla-extension"
334 : : #elif defined(STRICT_GNUC)
335 : : #pragma GCC diagnostic push
336 : : #pragma GCC diagnostic ignored "-Wvla"
337 : : #endif
338 : :
339 [ + + ]: 974720 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
340 : 969355 : const auto N = triinpoel.data() + e*3;
341 : :
342 : 969355 : auto rA = U(N[0],0);
343 : 969355 : auto ruA = U(N[0],1);
344 : 969355 : auto rvA = U(N[0],2);
345 : 969355 : auto rwA = U(N[0],3);
346 : 969355 : auto reA = U(N[0],4);
347 : :
348 : 969355 : auto rB = U(N[1],0);
349 : 969355 : auto ruB = U(N[1],1);
350 : 969355 : auto rvB = U(N[1],2);
351 : 969355 : auto rwB = U(N[1],3);
352 : 969355 : auto reB = U(N[1],4);
353 : :
354 : 969355 : auto rC = U(N[2],0);
355 : 969355 : auto ruC = U(N[2],1);
356 : 969355 : auto rvC = U(N[2],2);
357 : 969355 : auto rwC = U(N[2],3);
358 : 969355 : auto reC = U(N[2],4);
359 : :
360 : : const std::array< tk::real, 3 >
361 : 969355 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
362 [ + + ]: 969355 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
363 : : auto [nx,ny,nz] = tk::cross( ba, ca ); // 2A
364 : 969355 : nx /= 12.0;
365 : 969355 : ny /= 12.0;
366 : 969355 : nz /= 12.0;
367 : :
368 [ + + ]: 969355 : tk::real p, vn, f[ncomp][3];
369 : 969355 : const auto sym = besym.data() + e*3;
370 : :
371 [ + + ]: 969355 : p = eos::pressure( reA - 0.5*(ruA*ruA + rvA*rvA + rwA*rwA)/rA );
372 [ + + ]: 969355 : vn = sym[0] ? 0.0 : (nx*ruA + ny*rvA + nz*rwA)/rA;
373 : : // flow
374 : 969355 : f[0][0] = rA*vn;
375 : 969355 : f[1][0] = ruA*vn + p*nx;
376 : 969355 : f[2][0] = rvA*vn + p*ny;
377 : 969355 : f[3][0] = rwA*vn + p*nz;
378 : 969355 : f[4][0] = (reA + p)*vn;
379 : : // scalar
380 [ + + ]: 1113235 : for (std::size_t c=5; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
381 : :
382 [ + + ]: 969355 : p = eos::pressure( reB - 0.5*(ruB*ruB + rvB*rvB + rwB*rwB)/rB );
383 [ + + ]: 969355 : vn = sym[1] ? 0.0 : (nx*ruB + ny*rvB + nz*rwB)/rB;
384 : : // flow
385 : 969355 : f[0][1] = rB*vn;
386 : 969355 : f[1][1] = ruB*vn + p*nx;
387 : 969355 : f[2][1] = rvB*vn + p*ny;
388 : 969355 : f[3][1] = rwB*vn + p*nz;
389 : 969355 : f[4][1] = (reB + p)*vn;
390 : : // scalar
391 [ + + ]: 1113235 : for (std::size_t c=5; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
392 : :
393 [ + + ]: 969355 : p = eos::pressure( reC - 0.5*(ruC*ruC + rvC*rvC + rwC*rwC)/rC );
394 [ + + ]: 969355 : vn = sym[2] ? 0.0 : (nx*ruC + ny*rvC + nz*rwC)/rC;
395 : : // flow
396 : 969355 : f[0][2] = rC*vn;
397 : 969355 : f[1][2] = ruC*vn + p*nx;
398 : 969355 : f[2][2] = rvC*vn + p*ny;
399 : 969355 : f[3][2] = rwC*vn + p*nz;
400 : 969355 : f[4][2] = (reC + p)*vn;
401 : : // scalar
402 [ + + ]: 1113235 : for (std::size_t c=5; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
403 : :
404 [ + + ]: 5960010 : for (std::size_t c=0; c<ncomp; ++c) {
405 : 4990655 : auto fab = (f[c][0] + f[c][1])/4.0;
406 : 4990655 : auto fbc = (f[c][1] + f[c][2])/4.0;
407 : 4990655 : auto fca = (f[c][2] + f[c][0])/4.0;
408 : 4990655 : R(N[0],c) += fab + fca + f[c][0];
409 : 4990655 : R(N[1],c) += fab + fbc + f[c][1];
410 : 4990655 : R(N[2],c) += fbc + fca + f[c][2];
411 : : }
412 [ + + ]: 969355 : }
413 : :
414 : : #if defined(__clang__)
415 : : #pragma clang diagnostic pop
416 : : #elif defined(STRICT_GNUC)
417 : : #pragma GCC diagnostic pop
418 : : #endif
419 : 5365 : }
420 : :
421 : : void
422 : 5365 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
423 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
424 : : const std::array< std::vector< tk::real >, 3 >& coord,
425 : : const std::vector< std::size_t >& triinpoel,
426 : : const std::vector< std::uint8_t >& besym,
427 : : tk::real t,
428 : : tk::real dt,
429 : : const std::vector< tk::real >& tp,
430 : : const std::vector< tk::real >& dtp,
431 : : const tk::Fields& U,
432 : : tk::Fields& R )
433 : : // *****************************************************************************
434 : : // Compute right hand side
435 : : //! \param[in] dsupedge Domain superedges
436 : : //! \param[in] dsupint Domain superedge integrals
437 : : //! \param[in] coord Mesh node coordinates
438 : : //! \param[in] triinpoel Boundary face connectivity
439 : : //! \param[in] besym Boundary element symmetry BC flags
440 : : //! \param[in] t Physical time
441 : : //! \param[in] dt Physical time size
442 : : //! \param[in] tp Phisical time step size for each mesh node (if steady state)
443 : : //! \param[in] dtp Time step size for each mesh node (if steady state)
444 : : //! \param[in] U Unknowns/solution vector in mesh nodes
445 : : //! \param[in,out] R Right-hand side vector computed
446 : : // *****************************************************************************
447 : : {
448 : : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
449 : : "vector at recent time step incorrect" );
450 : : Assert( R.nunk() == coord[0].size(),
451 : : "Number of unknowns and/or number of components in right-hand "
452 : : "side vector incorrect" );
453 : :
454 : : R.fill( 0.0 );
455 : 5365 : advdom( dsupedge, dsupint, coord, t, dt, tp, dtp, U, R );
456 : 5365 : advbnd( triinpoel, coord, besym, U, R );
457 : 5365 : }
458 : :
459 : : } // zalesak::
|
