test_spaceoperator.cc

#include "function_library.hh"
#include "projection.hh"
#include "error.hh"
#include "mygrid.hh"
#include "finite_difference.hh"
#include "dofvector.hh"
#include "dataoutput.hh"
#include "operators/fv_hyperbolic.hh"

#include <iostream>
#include <sstream>

int main(int argc, const char *argv[])
{
  /*********************************************************************
   * instantiate all model objects which are independent of the grid,  *
   * i.e. the discretization                                           *
   ********************************************************************/
  typedef DefaultDofVectorTraits< double >                           DofVectorTraits;
  typedef DefaultDofVector< DofVectorTraits >                        DofVectorType;

  // instantiate model
  typedef TestDiscreteSpaceOperatorModel<TestFunc1> Testmodel;
  Testmodel model;
  // instantiate function and derivative from function library...
  TestFunc1      func;
  TestDerivFunc1 derivFunc;

  // instantiate the lax friedrichs flux
  const LaxFriedrichsFlux<Testmodel::FluxFunctionType> numFlux(model.f, 0);

  // initial number of grid cells
  int numintervals = 10;
  // maximum number of eoc steps
  unsigned int max_eocsteps;

  // read the maximum number of eocsteps from the command line...
  if (argc == 2)
  {
    std::istringstream iss(argv[1]);
    iss >> max_eocsteps;
  }
  else
  {
    // ... or set it to 5 by default
    max_eocsteps = 5;
  }


  // do the EOC loop
  for (unsigned int eocstep = 0; eocstep < max_eocsteps; eocstep++)
  {
    double oldl2error;
    double eoc;

    // instantiate grid on interval [0,3] with #numintervals grid cells
    Grid g(0, 3, numintervals);

    // instantiate discrete space operator
    SpaceOperator::FV::HyperbolicDefault<Testmodel> op(numFlux, model, g);

    // instantiate projection operator
    Projection projection(g);

    // project function, write result to dof vector arg
    DofVectorType arg(g);
    projection.apply(func, arg);
    DataOutput :: GnuplotOneD goned_arg(g, "function");
    goned_arg.write(arg, 0);

    // apply space operator, write result to dof vector res
    DofVectorType res(g);
    op.apply(arg, res);
    DataOutput :: GnuplotOneD goned_res(g, "res");
    goned_res.write(res, 0);

    // project derivative of function, write result to dof vector exact
    DofVectorType exact(g);
    projection.apply(derivFunc, exact);
    DataOutput :: GnuplotOneD goned_exact(g, "exact");
    goned_exact.write(exact, 0);

    // compute l2 error between exact and computed derivative...
    L2Error l2error(g);
    double l2err = l2error.apply(res, exact);
    // ... and print it out
    std::cout << "L2-Error: " << l2err << std::endl;

    // if we are in the second or later eoc step...
    if (eocstep > 0)
    {
      // ... we can compute an eoc and print it
      eoc = log(oldl2error/l2err)/log(2);
      std::cout << "EOC: " << eoc << std::endl;
    }
    oldl2error = l2err;

    // increase the number of grid cells for next run
    numintervals = numintervals * 2;
  }
  return 0;
}