1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
///
/// =========================================================================
//! @examplefile navier_stokes_dg1.icc The Navier-Stokes equations with the discontinuous Galerkin method -- class body
navier_stokes_dg::navier_stokes_dg (
Float Re1, const geo& omega, string approx)
: Re(Re1), Xh(), Qh(), iopt(), a0(), b(), c(), mu(), mp(), lh0(), lh(), kh(),
pmu(), pmp(), a1(), stokes1()
{
Xh = space (omega, approx, "vector");
Qh = space (omega, approx);
iopt.set_family(integrate_option::gauss);
iopt.set_order(2*Xh.degree()+1);
stokes_dirichlet_dg (Xh, Qh, a0, b, c, mp, lh0, kh, iopt);
trial u (Xh); test v (Xh);
lh = lh0 + Re*inertia_fix_rhs (v, iopt);
mu = integrate (dot(u,v), iopt);
pmu = problem (mu);
pmp = problem (mp);
}
navier_stokes_dg::value_type
navier_stokes_dg::initial (string restart) const {
value_type xh = { field(Xh, 0), field(Qh, 0) };
Float Re0 = 0;
if (restart == "") {
problem_mixed stokes0 (a0, b, c);
stokes0.set_metric (mp);
stokes0.solve (lh0, kh, xh[0], xh[1]);
} else {
idiststream in (restart);
in >> catchmark("Re") >> Re0
>> catchmark("u") >> xh[0]
>> catchmark("p") >> xh[1];
check_macro (xh[1].get_space() == Qh, "unexpected "
<< xh[0].get_space().name() << " approximation in file \""
<< restart << "\" (" << Xh.name() << " expected)");
}
derr << "# continuation: from Re=" << Re0 << " to " << Re << endl;
return xh;
}
navier_stokes_dg::value_type
navier_stokes_dg::residue (const value_type& xh) const {
trial u (Xh); test v (Xh);
form a = a0 + Re*inertia(xh[0], u, v, iopt);
value_type mrh = { a*xh[0] + b.trans_mult(xh[1]) - lh,
b*xh[0] - c*xh[1] - kh};
return mrh;
}
void navier_stokes_dg::update_derivative (const value_type& xh) const {
trial u (Xh); test v (Xh);
a1 = a0 + Re*(inertia(xh[0], u, v, iopt) + inertia(u, xh[0], v, iopt));
stokes1 = problem_mixed (a1, b, c);
stokes1.set_metric (mp);
}
navier_stokes_dg::value_type
navier_stokes_dg::derivative_solve (const value_type& mrh) const {
value_type delta_xh = { field(Xh, 0), field(Qh, 0) };
stokes1.solve (mrh[0], mrh[1], delta_xh[0], delta_xh[1]);
return delta_xh;
}
navier_stokes_dg::value_type
navier_stokes_dg::derivative_trans_mult (const value_type& mrh) const {
value_type rh = { field (Xh), field (Qh) };
pmu.solve (mrh[0], rh[0]);
pmp.solve (mrh[1], rh[1]);
value_type mgh = { a1.trans_mult(rh[0]) + b.trans_mult(rh[1]),
b*rh[0] - c*rh[1] };
return mgh;
}
|
