limiter.cc
Go to the documentation of this file.
// minmod_TVB limiter for hyperbolic nonlinear problems
// approximated by discontinuous Galerkin methods
//
#include "rheolef/limiter.h"
#include "rheolef/integrate.h"
namespace rheolef {
// ----------------
// minmod functions
// ----------------
namespace details {
template <class T>
inline
T
minmod (const T& a, const T& b) {
// avoid a=1e-17 and b=-1 => minmod(a,b)=0 instead of b
static T epsilon = std::numeric_limits<T>::epsilon();
static T tol = sqrt(epsilon);
if (fabs(a) < tol) return 0;
if (fabs(b) < tol) return 0;
T sgn_a = (a >= 0) ? 1 : -1;
T sgn_b = (b >= 0) ? 1 : -1;
T res = (sgn_a == sgn_b) ? sgn_a*min(fabs(a),fabs(b)) : 0;
trace_macro ("minmod("<<a<<","<<b<<") = " << res);
return res;
}
template <class T>
inline
T
minmod_tvb (const T& yield_a, const T& a, const T& b) {
T res = (fabs(a) <= yield_a) ? a : minmod (a,b);
trace_macro ("minmod_tdb("<<yield_a<<","<<a<<","<<b<<") = " << res);
return res;
}
} // namespace details
// ----------------
// limiter function
// ----------------
template <class T, class M>
field_basic<T,M>
limiter (
const field_basic<T,M>& uh,
const T& bar_g_S, // TODO: general boundary condition
const limiter_option& opt)
{
if (! opt.active) return uh;
typedef typename field_basic<T,M>::size_type size_type;
T epsilon = std::numeric_limits<T>::epsilon();
T tol = sqrt(epsilon);
const geo_basic<T,M>& omega = uh.get_geo();
omega.neighbour_guard();
check_macro (uh.get_space().get_basis().is_discontinuous(),
"argument may be a discontinuous approximation: "
<< uh.get_space().name() << " founded");
size_type k = uh.get_space().degree();
if (k == 0) return uh;
check_macro (k == 1, "order = " << k << " not yet supported");
size_type map_d = omega.map_dimension();
check_macro (map_d == 1, "dimension map_d = " << map_d << " not yet supported");
uh.dis_dof_update();
// 0. measure of each elements
space_basic<T,M> X0h (omega, "P0");
test_basic<T,M,details::vf_tag_01> v (X0h);
field_basic<T,M> meas_K_var = integrate (v);
meas_K_var.dis_dof_update();
const field_basic<T,M>& meas_K = meas_K_var;
constexpr size_type ns_loc_max = 8; // subgeo (hexa -> 8 quadri)
constexpr size_type nss_loc_max = 4; // subsubgeo (quadri -> 4 edges)
std::array<bool,ns_loc_max> is_on_boundary;
std::array<bool,ns_loc_max> is_upstream;
std::array<point_basic<T>,ns_loc_max> xK1;
size_type index [ns_loc_max][nss_loc_max];
T alpha [ns_loc_max][nss_loc_max];
std::array<T,ns_loc_max> tilde, delta, Delta, hat_Delta;
std::array<geo_element::orientation_type,ns_loc_max> orient;
field_basic<T,M> vh (uh.get_space(), 0);
for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
const geo_element& K = omega.get_geo_element (map_d, ie);
trace_macro ("K"<<K.dis_ie()<<" ios_dis_ie="<<K.ios_dis_ie()<<"...");
// --------------
// geometric data
// --------------
// 1. compute the barycenter of K
trace_macro ("K"<<K.dis_ie()<<": barycenter...");
size_type nv_loc = K.size();
point_basic<T> xK;
for (size_type iv_loc = 0; iv_loc < nv_loc; ++iv_loc) {
size_type dis_iv = K[iv_loc];
xK += omega.dis_node (dis_iv);
}
xK /= T(int(nv_loc));
size_type ndof_per_K = 2; // TODO: 1D only
trace_macro ("K"<<K.dis_ie()<<": index and alpha...");
// 2. compute the geometric data
for (size_type is_loc = 0, ns_loc = K.n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
trace_macro ("K"<<K.dis_ie()<<": is_loc="<<is_loc<<"...");
size_type dis_is = (map_d == 1) ? K[is_loc] : ((map_d == 2) ? K.edge(is_loc) : K.face(is_loc));
trace_macro ("K"<<K.dis_ie()<<": dis_is="<<dis_is);
const geo_element& S = omega.dis_get_geo_element (map_d-1, dis_is);
trace_macro ("K"<<K.dis_ie()<<": S.dis_ie="<<S.dis_ie());
trace_macro ("K"<<K.dis_ie()<<": S.master0="<<S.master(0));
trace_macro ("K"<<K.dis_ie()<<": S.master1="<<S.master(1));
// get Ki = neighbour of K across Si:
size_type dis_ie1 = (S.master(0) != K.dis_ie()) ? S.master(0) : S.master(1);
trace_macro ("K"<<K.dis_ie()<<": dis_ie1="<<dis_ie1);
is_on_boundary [is_loc] = (dis_ie1 == std::numeric_limits<size_type>::max());
if (is_on_boundary[is_loc]) {
index [is_loc][0] = is_loc; // TODO: 1D only -> test determinant sign
alpha [is_loc][0] = 1;
continue;
}
const geo_element& K1 = omega.dis_get_geo_element (map_d, dis_ie1);
trace_macro ("K"<<K.dis_ie()<<": K1.dis_ie="<<K1.dis_ie());
// 2.1. compute the barycenter of Ki
xK1 [is_loc] = {0,0,0};
size_type nv1_loc = K1.size();
for (size_type iv1_loc = 0; iv1_loc < nv1_loc; ++iv1_loc) {
size_type dis_iv1 = K1[iv1_loc];
trace_macro ("K"<<K.dis_ie()<<": dis_iv1="<<dis_iv1);
xK1 [is_loc] += omega.dis_node (dis_iv1);
}
xK1 [is_loc] /= T(int(nv1_loc));
// 2.2. compute index(i,k) and alpha(i,k)
trace_macro ("K"<<K.dis_ie()<<": index and alpha...");
index [is_loc][0] = is_loc; // TODO: 1D only -> test determinant sign
alpha [is_loc][0] = meas_K.dof(ie)/(meas_K.dis_dof(dis_ie1) + meas_K.dof(ie));
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: alpha="<<alpha[is_loc][0]);
}
// -------------------------
// uh-dependent computations
// -------------------------
// 1. compute the average value of u on K
trace_macro ("K"<<K.dis_ie()<<": bar_u_K...");
T bar_u_K = 0;
for (size_type iv_loc = 0; iv_loc < nv_loc; ++iv_loc) {
size_type idof = ndof_per_K*ie + iv_loc;
bar_u_K += uh.dof (idof);
}
bar_u_K /= nv_loc;
trace_macro ("K"<<K.dis_ie()<<": bar_u_K="<<bar_u_K);
T hK = pow (meas_K.dof(ie), 1./map_d); // used only by minmod_tvb()
// 2. tilde, delta, Delta and hat_Delta
T sum_Delta_plus = 0,
sum_Delta_minus = 0;
for (size_type is_loc = 0, ns_loc = K.n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
// get Ki = neighbour of K across Si:
size_type dis_is = (map_d == 1) ? K[is_loc] : ((map_d == 2) ? K.edge(is_loc) : K.face(is_loc));
const geo_element& S = omega.dis_get_geo_element (map_d-1, dis_is);
orient[is_loc] = 1;
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: orient="<<orient[is_loc]);
size_type dis_ie1 = (S.master(0) != K.dis_ie()) ? S.master(0) : S.master(1);
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: dis_ie1="<<dis_ie1);
is_on_boundary [is_loc] = (dis_ie1 == std::numeric_limits<size_type>::max());
is_upstream[is_loc] = is_on_boundary [is_loc] && (K.ios_dis_ie() == 0); // TODO: 2D-only and done with mkgeo_grid...
// 2.1. tilde
T bar_u_S;
if (false && is_upstream [is_loc]) {
bar_u_S = bar_g_S;
} else {
size_type idof = ndof_per_K*ie + is_loc; // 1D-only
bar_u_S = uh.dof (idof); // TODO: 1D ony -> average over S vertices
}
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: bar_u_S="<<bar_u_S);
tilde [is_loc] = orient[is_loc]*(bar_u_S - bar_u_K);
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: tilde="<<tilde[is_loc]);
// 2.2. delta
T bar_u_K1 = 0;
if (is_on_boundary [is_loc]) {
if (is_upstream [is_loc]) {
bar_u_K1 = bar_g_S;
} else {
bar_u_K1 = bar_u_K;
}
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: bar_u_K1="<<bar_u_K1);
} else {
const geo_element& K1 = omega.dis_get_geo_element (map_d, dis_ie1);
size_type nv1_loc = K1.size();
for (size_type iv1_loc = 0; iv1_loc < nv1_loc; ++iv1_loc) {
size_type dis_iv1 = K1[iv1_loc];
size_type dis_idof = ndof_per_K*K1.dis_ie() + iv1_loc;
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: uh on K1...");
bar_u_K1 += uh.dis_dof (dis_idof);
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: uh on K1 done");
}
bar_u_K1 /= nv1_loc;
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: bar_u_K1="<<bar_u_K1);
}
delta [is_loc] = orient[is_loc]*alpha[is_loc][0]*(bar_u_K1 - bar_u_K); // TODO: sum_k
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: delta="<<delta[is_loc]);
// 2.3. Delta
Delta [is_loc] = details::minmod_tvb (opt.M*sqr(hK), tilde[is_loc], opt.theta*delta[is_loc]);
sum_Delta_plus += max (T(0.), Delta [is_loc]);
sum_Delta_minus += max (T(0.), -Delta [is_loc]);
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: Delta="<<Delta[is_loc]);
}
// 2.4. hat_Delta
T r = (fabs(sum_Delta_plus) > tol) ? sum_Delta_minus/sum_Delta_plus : 0;
trace_macro ("K"<<K.dis_ie()<<" sum_Delta_plus ="<<sum_Delta_plus);
trace_macro ("K"<<K.dis_ie()<<" sum_Delta_minus="<<sum_Delta_minus);
trace_macro ("K"<<K.dis_ie()<<" r ="<<r);
for (size_type is_loc = 0, ns_loc = K.n_subgeo(map_d-1); is_loc < ns_loc; ++is_loc) {
hat_Delta [is_loc] = (r == 0) ? Delta [is_loc]
: min(T(1.), r)*max(T(0.), Delta [is_loc])
- min(T(1.),1./r)*max(T(0.), -Delta [is_loc]);
trace_macro ("K"<<K.dis_ie()<<"["<<is_loc<<"]: hat_Delta="<<hat_Delta[is_loc]);
}
// 2.5. vh = limiter(uh)
for (size_type iv_loc = 0, nv_loc = K.size(); iv_loc < nv_loc; ++iv_loc) {
size_type idof = ndof_per_K*ie + iv_loc; // TODO: P1d-only idof -> interpolate P1d at others Pkd local_xdofs
size_type is_loc = iv_loc; // TODO: P1d-only
if (false && is_on_boundary [is_loc] && ! is_upstream[is_loc]) {
vh.dof (idof) = uh.dof (idof);
} else {
vh.dof (idof) = bar_u_K + orient[is_loc]*hat_Delta [is_loc]; // TODO: 1D-only: interpolate at side-based Lagrange basis phi_is
}
trace_macro ("K"<<K.dis_ie()<<"["<<iv_loc<<"]: idof="<<idof<<", x="<<ptos(uh.get_space().xdof(idof),1)
<< ", u="<<uh.dof (idof)<<" -> " << vh.dof (idof));
}
}
vh.dis_dof_update();
return vh;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M) \
template field_basic<T,M> limiter ( \
const field_basic<T,M>&, \
const T&, \
const limiter_option&);
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
rheolef::field_basic< T, M >
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::field_basic::dis_dof
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:376
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::field_basic::dis_dof_update
void dis_dof_update() const
Definition: field.cc:410
rheolef::field_basic::dof
T & dof(size_type idof)
Definition: field.h:658
rheolef::geo_basic< T, M >
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::edge
size_type edge(size_type i) const
Definition: geo_element.h:209
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element::ios_dis_ie
size_type ios_dis_ie() const
Definition: geo_element.h:164
rheolef::geo_element::master
size_type master(bool i) const
Definition: geo_element.h:165
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::point_basic
Definition: point.h:87
rheolef::space_basic< T, M >
rheolef::test_basic< T, M, details::vf_tag_01 >
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::distributed
distributed
Definition: asr.cc:228
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
delta
Float delta(Float f, Float g)
Definition: mosolov_error_yield_surface.cc:29
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::details::minmod_tvb
T minmod_tvb(const T &yield_a, const T &a, const T &b)
Definition: limiter.cc:51
rheolef::details::minmod
T minmod(const T &a, const T &b)
Definition: limiter.cc:36
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
This file is part of Rheolef.
Definition: compiler_eigen.h:37
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::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::ptos
std::string ptos(const point_basic< T > &x, int d=3)
Definition: point.h:414
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
a
Definition: diffusion_isotropic.h:25
rheolef::limiter_option
see the limiter page for the full documentation
Definition: limiter.h:72
rheolef::limiter_option::active
bool active
Definition: limiter.h:73
rheolef::limiter_option::M
Float M
Definition: limiter.h:75
rheolef::limiter_option::theta
Float theta
Definition: limiter.h:74
epsilon
Float epsilon
Definition: transmission_error.cc:25