Rheolef  7.1
an efficient C++ finite element environment
limiter.cc
Go to the documentation of this file.
1 // 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 
27 namespace rheolef {
28 
29 // ----------------
30 // minmod functions
31 // ----------------
32 namespace details {
33 template <class T>
34 inline
35 T
36 minmod (const T& a, const T& b) {
37  // avoid a=1e-17 and b=-1 => minmod(a,b)=0 instead of b
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 }
48 template <class T>
49 inline
50 T
51 minmod_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 // ----------------
63 template <class T, class M>
64 field_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;
71  typedef typename field_basic<T,M>::size_type size_type;
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();
112  point_basic<T> xK;
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) \
253 template field_basic<T,M> limiter ( \
254  const field_basic<T,M>&, \
255  const T&, \
256  const limiter_option&);
257 
258 _RHEOLEF_instanciation(Float,sequential)
259 #ifdef _RHEOLEF_HAVE_MPI
261 #endif // _RHEOLEF_HAVE_MPI
262 
263 } // namespace rheolef
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::point_basic
Definition: point.h:87
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::geo_element::master
size_type master(bool i) const
Definition: geo_element.h:165
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::integrate
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:202
rheolef::limiter_option
see the limiter page for the full documentation
Definition: limiter.h:72
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef::field_basic::dis_dof
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:376
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::test_basic< T, M, details::vf_tag_01 >
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::details::minmod_tvb
T minmod_tvb(const T &yield_a, const T &a, const T &b)
Definition: limiter.cc:51
a
Definition: diffusion_isotropic.h:25
rheolef::field_basic
Definition: field.h:235
rheolef::field_basic::dof
T & dof(size_type idof)
Definition: field.h:658
rheolef::limiter
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
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::limiter_option::active
bool active
Definition: limiter.h:73
rheolef::details::minmod
T minmod(const T &a, const T &b)
Definition: limiter.cc:36
Float
see the Float page for the full documentation
rheolef::limiter_option::M
Float M
Definition: limiter.h:75
rheolef::ptos
std::string ptos(const point_basic< T > &x, int d=3)
Definition: point.h:414
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
epsilon
Float epsilon
Definition: transmission_error.cc:25
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209
delta
Float delta(Float f, Float g)
Definition: mosolov_error_yield_surface.cc:29
rheolef::field_basic::size_type
std::size_t size_type
Definition: field.h:239
rheolef::geo_element::ios_dis_ie
size_type ios_dis_ie() const
Definition: geo_element.h:164
T
Expr1::float_type T
Definition: field_expr.h:261
rheolef::limiter_option::theta
Float theta
Definition: limiter.h:74
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111