convect_error.cc

Convection-diffusion equation by the method of characteristics – error analysis.

Convection-diffusion equation by the method of characteristics – error analysis

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "rotating-hill.h"
int main (int argc, char **argv) {
  environment rheolef (argc,argv);
  Float tol = (argc > 1) ? atof(argv[1]) : 1e-10;
  Float nu;
  din >> catchmark("nu") >> nu;
  branch get ("t","phi");
  branch put ("t","phi_h","pi_h_phi");
  derr << "# t\terror_l2\terror_linf" << endl;
  field phi_h;
  Float err_l2_l2 = 0;
  Float err_linf_linf = 0;
  for (Float t = 0, t_prec = 0; din >> get (t, phi_h); t_prec = t) {
    const space& Xh = phi_h.get_space();
    size_t d = Xh.get_geo().dimension();
    field pi_h_phi = lazy_interpolate (Xh, phi(d,nu,t));
    trial phi (Xh); test psi (Xh);
    form m = integrate (phi*psi);
    field eh = phi_h - pi_h_phi;
    Float err_l2 = sqrt(m(eh,eh));
    Float err_linf = eh.max_abs();
    err_l2_l2 += sqr(err_l2)*(t - t_prec);
    err_linf_linf = max(err_linf_linf, err_linf);
    dout << put (t, phi_h, pi_h_phi);
    derr << t << "\t" << err_l2 << "\t" << err_linf << endl;
  }
  derr << "# error_l2_l2     = " << sqrt(err_l2_l2) << endl;
  derr << "# error_linf_linf = " << err_linf_linf << endl;
  return (err_linf_linf <= tol) ? 0 : 1;
}