Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
limiter.cc
Go to the documentation of this file.
1
21// minmod_TVB limiter for hyperbolic nonlinear problems
22// approximated by discontinuous Galerkin methods
23//
24#include "rheolef/limiter.h"
25#include "rheolef/integrate.h"
26
27namespace rheolef {
28
29// ----------------
30// minmod functions
31// ----------------
32namespace details {
33template <class T>
34inline
35T
36minmod (const T& a, const T& b) {
37 // avoid a=1e-17 and b=-1 => minmod(a,b)=0 instead of b
38 static T epsilon = std::numeric_limits<T>::epsilon();
39 static T tol = sqrt(epsilon);
40 if (fabs(a) < tol) return 0;
41 if (fabs(b) < tol) return 0;
42 T sgn_a = (a >= 0) ? 1 : -1;
43 T sgn_b = (b >= 0) ? 1 : -1;
44 T res = (sgn_a == sgn_b) ? sgn_a*min(fabs(a),fabs(b)) : 0;
45 trace_macro ("minmod("<<a<<","<<b<<") = " << res);
46 return res;
47}
48template <class T>
49inline
50T
51minmod_tvb (const T& yield_a, const T& a, const T& b) {
52 T res = (fabs(a) <= yield_a) ? a : minmod (a,b);
53 trace_macro ("minmod_tdb("<<yield_a<<","<<a<<","<<b<<") = " << res);
54 return res;
55}
56
57} // namespace details
58
59// ----------------
60// limiter function
61// ----------------
63template <class T, class M>
64field_basic<T,M>
66 const field_basic<T,M>& uh,
67 const T& bar_g_S, // TODO: general boundary condition
68 const limiter_option& opt)
69{
70 if (! opt.active) return uh;
72 T epsilon = std::numeric_limits<T>::epsilon();
73 T tol = sqrt(epsilon);
74 const geo_basic<T,M>& omega = uh.get_geo();
75 omega.neighbour_guard();
76 check_macro (uh.get_space().get_basis().is_discontinuous(),
77 "argument may be a discontinuous approximation: "
78 << uh.get_space().name() << " founded");
79 size_type k = uh.get_space().degree();
80 if (k == 0) return uh;
81 check_macro (k == 1, "order = " << k << " not yet supported");
82 size_type map_d = omega.map_dimension();
83 check_macro (map_d == 1, "dimension map_d = " << map_d << " not yet supported");
84 uh.dis_dof_update();
85
86 // 0. measure of each elements
87 space_basic<T,M> X0h (omega, "P0");
89 field_basic<T,M> meas_K_var = integrate (v);
90 meas_K_var.dis_dof_update();
91 const field_basic<T,M>& meas_K = meas_K_var;
92
93 constexpr size_type ns_loc_max = 8; // subgeo (hexa -> 8 quadri)
94 constexpr size_type nss_loc_max = 4; // subsubgeo (quadri -> 4 edges)
95 std::array<bool,ns_loc_max> is_on_boundary;
96 std::array<bool,ns_loc_max> is_upstream;
97 std::array<point_basic<T>,ns_loc_max> xK1;
98 size_type index [ns_loc_max][nss_loc_max];
99 T alpha [ns_loc_max][nss_loc_max];
100 std::array<T,ns_loc_max> tilde, delta, Delta, hat_Delta;
101 std::array<geo_element::orientation_type,ns_loc_max> orient;
102 field_basic<T,M> vh (uh.get_space(), 0);
103 for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
104 const geo_element& K = omega.get_geo_element (map_d, ie);
105 trace_macro ("K"<<K.dis_ie()<<" ios_dis_ie="<<K.ios_dis_ie()<<"...");
106 // --------------
107 // geometric data
108 // --------------
109 // 1. compute the barycenter of K
110 trace_macro ("K"<<K.dis_ie()<<": barycenter...");
111 size_type nv_loc = K.size();
113 for (size_type iv_loc = 0; iv_loc < nv_loc; ++iv_loc) {
114 size_type dis_iv = K[iv_loc];
115 xK += omega.dis_node (dis_iv);
116 }
117 xK /= T(int(nv_loc));
118 size_type ndof_per_K = 2; // TODO: 1D only
119 trace_macro ("K"<<K.dis_ie()<<": index and alpha...");
120 // 2. compute the geometric data
121 for (size_type is_loc = 0, ns_loc = K.n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
122 trace_macro ("K"<<K.dis_ie()<<": is_loc="<<is_loc<<"...");
123 size_type dis_is = (map_d == 1) ? K[is_loc] : ((map_d == 2) ? K.edge(is_loc) : K.face(is_loc));
124 trace_macro ("K"<<K.dis_ie()<<": dis_is="<<dis_is);
125 const geo_element& S = omega.dis_get_geo_element (map_d-1, dis_is);
126 trace_macro ("K"<<K.dis_ie()<<": S.dis_ie="<<S.dis_ie());
127 trace_macro ("K"<<K.dis_ie()<<": S.master0="<<S.master(0));
128 trace_macro ("K"<<K.dis_ie()<<": S.master1="<<S.master(1));
129 // get Ki = neighbour of K across Si:
130 size_type dis_ie1 = (S.master(0) != K.dis_ie()) ? S.master(0) : S.master(1);
131 trace_macro ("K"<<K.dis_ie()<<": dis_ie1="<<dis_ie1);
132 is_on_boundary [is_loc] = (dis_ie1 == std::numeric_limits<size_type>::max());
133 if (is_on_boundary[is_loc]) {
134 index [is_loc][0] = is_loc; // TODO: 1D only -> test determinant sign
135 alpha [is_loc][0] = 1;
136 continue;
137 }
138 const geo_element& K1 = omega.dis_get_geo_element (map_d, dis_ie1);
139 trace_macro ("K"<<K.dis_ie()<<": K1.dis_ie="<<K1.dis_ie());
140 // 2.1. compute the barycenter of Ki
141 xK1 [is_loc] = {0,0,0};
142 size_type nv1_loc = K1.size();
143 for (size_type iv1_loc = 0; iv1_loc < nv1_loc; ++iv1_loc) {
144 size_type dis_iv1 = K1[iv1_loc];
145 trace_macro ("K"<<K.dis_ie()<<": dis_iv1="<<dis_iv1);
146 xK1 [is_loc] += omega.dis_node (dis_iv1);
147 }
148 xK1 [is_loc] /= T(int(nv1_loc));
149 // 2.2. compute index(i,k) and alpha(i,k)
150 trace_macro ("K"<<K.dis_ie()<<": index and alpha...");
151 index [is_loc][0] = is_loc; // TODO: 1D only -> test determinant sign
152 alpha [is_loc][0] = meas_K.dof(ie)/(meas_K.dis_dof(dis_ie1) + meas_K.dof(ie));
153 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: alpha="<<alpha[is_loc][0]);
154 }
155 // -------------------------
156 // uh-dependent computations
157 // -------------------------
158 // 1. compute the average value of u on K
159 trace_macro ("K"<<K.dis_ie()<<": bar_u_K...");
160 T bar_u_K = 0;
161 for (size_type iv_loc = 0; iv_loc < nv_loc; ++iv_loc) {
162 size_type idof = ndof_per_K*ie + iv_loc;
163 bar_u_K += uh.dof (idof);
164 }
165 bar_u_K /= nv_loc;
166 trace_macro ("K"<<K.dis_ie()<<": bar_u_K="<<bar_u_K);
167 T hK = pow (meas_K.dof(ie), 1./map_d); // used only by minmod_tvb()
168 // 2. tilde, delta, Delta and hat_Delta
169 T sum_Delta_plus = 0,
170 sum_Delta_minus = 0;
171 for (size_type is_loc = 0, ns_loc = K.n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
172 // get Ki = neighbour of K across Si:
173 size_type dis_is = (map_d == 1) ? K[is_loc] : ((map_d == 2) ? K.edge(is_loc) : K.face(is_loc));
174 const geo_element& S = omega.dis_get_geo_element (map_d-1, dis_is);
175 orient[is_loc] = 1;
176 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: orient="<<orient[is_loc]);
177 size_type dis_ie1 = (S.master(0) != K.dis_ie()) ? S.master(0) : S.master(1);
178 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: dis_ie1="<<dis_ie1);
179 is_on_boundary [is_loc] = (dis_ie1 == std::numeric_limits<size_type>::max());
180 is_upstream[is_loc] = is_on_boundary [is_loc] && (K.ios_dis_ie() == 0); // TODO: 2D-only and done with mkgeo_grid...
181 // 2.1. tilde
182 T bar_u_S;
183 if (false && is_upstream [is_loc]) {
184 bar_u_S = bar_g_S;
185 } else {
186 size_type idof = ndof_per_K*ie + is_loc; // 1D-only
187 bar_u_S = uh.dof (idof); // TODO: 1D ony -> average over S vertices
188 }
189 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: bar_u_S="<<bar_u_S);
190 tilde [is_loc] = orient[is_loc]*(bar_u_S - bar_u_K);
191 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: tilde="<<tilde[is_loc]);
192 // 2.2. delta
193 T bar_u_K1 = 0;
194 if (is_on_boundary [is_loc]) {
195 if (is_upstream [is_loc]) {
196 bar_u_K1 = bar_g_S;
197 } else {
198 bar_u_K1 = bar_u_K;
199 }
200 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: bar_u_K1="<<bar_u_K1);
201 } else {
202 const geo_element& K1 = omega.dis_get_geo_element (map_d, dis_ie1);
203 size_type nv1_loc = K1.size();
204 for (size_type iv1_loc = 0; iv1_loc < nv1_loc; ++iv1_loc) {
205 size_type dis_iv1 = K1[iv1_loc];
206 size_type dis_idof = ndof_per_K*K1.dis_ie() + iv1_loc;
207 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: uh on K1...");
208 bar_u_K1 += uh.dis_dof (dis_idof);
209 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: uh on K1 done");
210 }
211 bar_u_K1 /= nv1_loc;
212 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: bar_u_K1="<<bar_u_K1);
213 }
214 delta [is_loc] = orient[is_loc]*alpha[is_loc][0]*(bar_u_K1 - bar_u_K); // TODO: sum_k
215 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: delta="<<delta[is_loc]);
216 // 2.3. Delta
217 Delta [is_loc] = details::minmod_tvb (opt.M*sqr(hK), tilde[is_loc], opt.theta*delta[is_loc]);
218 sum_Delta_plus += max (T(0.), Delta [is_loc]);
219 sum_Delta_minus += max (T(0.), -Delta [is_loc]);
220 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: Delta="<<Delta[is_loc]);
221 }
222 // 2.4. hat_Delta
223 T r = (fabs(sum_Delta_plus) > tol) ? sum_Delta_minus/sum_Delta_plus : 0;
224 trace_macro ("K"<<K.dis_ie()<<" sum_Delta_plus ="<<sum_Delta_plus);
225 trace_macro ("K"<<K.dis_ie()<<" sum_Delta_minus="<<sum_Delta_minus);
226 trace_macro ("K"<<K.dis_ie()<<" r ="<<r);
227 for (size_type is_loc = 0, ns_loc = K.n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
228 hat_Delta [is_loc] = (r == 0) ? Delta [is_loc]
229 : min(T(1.), r)*max(T(0.), Delta [is_loc])
230 - min(T(1.),1./r)*max(T(0.), -Delta [is_loc]);
231 trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: hat_Delta="<<hat_Delta[is_loc]);
232 }
233 // 2.5. vh = limiter(uh)
234 for (size_type iv_loc = 0, nv_loc = K.size(); iv_loc < nv_loc; ++iv_loc) {
235 size_type idof = ndof_per_K*ie + iv_loc; // TODO: P1d-only idof -> interpolate P1d at others Pkd local_xdofs
236 size_type is_loc = iv_loc; // TODO: P1d-only
237 if (false && is_on_boundary [is_loc] && ! is_upstream[is_loc]) {
238 vh.dof (idof) = uh.dof (idof);
239 } else {
240 vh.dof (idof) = bar_u_K + orient[is_loc]*hat_Delta [is_loc]; // TODO: 1D-only: interpolate at side-based Lagrange basis phi_is
241 }
242 trace_macro ("K"<<K.dis_ie()<<"["<<iv_loc<<"]: idof="<<idof<<", x="<<ptos(uh.get_space().xdof(idof),1)
243 << ", u="<<uh.dof (idof)<<" -> " << vh.dof (idof));
244 }
245 }
246 vh.dis_dof_update();
247 return vh;
248}
249// ----------------------------------------------------------------------------
250// instanciation in library
251// ----------------------------------------------------------------------------
252#define _RHEOLEF_instanciation(T,M) \
253template field_basic<T,M> limiter ( \
254 const field_basic<T,M>&, \
255 const T&, \
256 const limiter_option&);
257
259#ifdef _RHEOLEF_HAVE_MPI
260_RHEOLEF_instanciation(Float,distributed)
261#endif // _RHEOLEF_HAVE_MPI
262
263} // namespace rheolef
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
const T & dis_dof(size_type dis_idof) const
Definition field.cc:377
T & dof(size_type idof)
Definition field.h:738
void dis_dof_update(const SetOp &=SetOp()) const
Definition field.h:763
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
std::size_t size_type
Definition field.h:225
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
size_type size() const
size_type ios_dis_ie() const
size_type master(bool i) const
size_type n_subgeo(size_type subgeo_dim) const
size_type dis_ie() const
the finite element space
Definition space.h:382
#define trace_macro(message)
Definition dis_macros.h:111
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
Float delta(Float f, Float g)
T minmod_tvb(const T &yield_a, const T &a, const T &b)
Definition limiter.cc:51
T minmod(const T &a, const T &b)
Definition limiter.cc:36
This file is part of Rheolef.
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&!is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result dummy=Result())
see the integrate page for the full documentation
Definition integrate.h:211
field_basic< T, M > limiter(const field_basic< T, M > &uh, const T &bar_g_S, const limiter_option &opt)
see the limiter page for the full documentation
Definition limiter.cc:65
std::string ptos(const point_basic< T > &x, int d=3)
Definition point.h:413
see the limiter page for the full documentation
Definition limiter.h:72
Float epsilon