zalesak_dg_error.cc

The Zalesak slotted disk benchmark – error analysis.

The Zalesak slotted disk benchmark – error analysis

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "zalesak.h"
Float heaviside (const Float& x) { return (x <= 0) ? 0 : 1; }
int main(int argc, char**argv) {
  environment rheolef (argc, argv);
  size_t subdivide = (argc > 1) ? atoi(argv[1]) : 1;
  Float tol        = (argc > 2) ? atof(argv[2]) : 1e-12;
  Float meas_gamma = phi_exact(0).perimeter();
  branch event ("t","phi");
  integrate_option iopt;
  iopt.set_family ("equispaced");
  iopt.set_order  (subdivide);
  dout << setprecision(numeric_limits<Float>::digits10)
       << "# meas_gamma = " << meas_gamma << endl
       << "# equispaced("<<subdivide<<")" << endl
       << "# t err_l1 v=(m-m0)/m0 m" << endl;
  bool first = true;
  Float t = 0, err_linf_l1 = 0, m0 = 0, v_max = 0;
  field phi_h;
  geo omega;
  while (din >> event (t, phi_h)) {  
    omega = phi_h.get_geo();
    Float err_l1 = integrate (omega, abs(compose(heaviside,phi_h)
                                       - compose(heaviside,phi_exact(t))), iopt);
    Float m = integrate (omega, 1. - compose(heaviside,phi_h), iopt);
    if (first) {
      m0 = m;
      first = false;
    }
    Float v = (m-m0)/m0;
    dout << t << " " << err_l1/meas_gamma 
              << " " << v << " " << m << endl;
    err_linf_l1 = max (err_linf_l1, err_l1/meas_gamma);
    v_max = max (v_max, abs(v));
  }
  dout << "# omega.size  = " << omega.size() << endl
       << "# omega.hmin  = " << omega.hmin() << endl
       << "# err_linf_l1 = " << err_linf_l1 << endl
       << "# v_max       = " << v_max << endl;
  return (err_linf_l1 < tol && v_max < tol) ? 0 : 1;
}