finitevolume.cc 3.36 KB
Newer Older
1 2
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
3
#include <config.h>               // know what grids are present
4 5 6
#include <iostream>               // for input/output to shell
#include <fstream>                // for input/output to files
#include <vector>                 // STL vector class
7
#include <dune/grid/common/mcmgmapper.hh> // mapper class
8
#include <dune/common/parallel/mpihelper.hh> // include mpi helper class
Robert Klöfkorn's avatar
Robert Klöfkorn committed
9 10

#include "vtkout.hh"
11
#include "transportproblem2.hh"
12 13 14 15 16 17 18 19 20 21 22
#include "initialize.hh"
#include "evolve.hh"

//===============================================================
// the time loop function working for all types of grids
//===============================================================

template<class G>
void timeloop (const G& grid, double tend)
{
  // make a mapper for codim 0 entities in the leaf grid
23 24
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G>
  mapper(grid, Dune::mcmgElementLayout());
25 26 27 28 29

  // allocate a vector for the concentration
  std::vector<double> c(mapper.size());

  // initialize concentration with initial values
30
  initialize(grid,mapper,c);                           /*@\label{fvc:init}@*/
31
  vtkout(grid,c,"concentration",0,0.0);
32 33 34 35

  // now do the time steps
  double t=0,dt;
  int k=0;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
36 37
  const double saveInterval = 0.1;
  double saveStep = 0.1;
38
  int counter = 1;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
39

40
  while (t<tend)                                       /*@\label{fvc:loop0}@*/
41
  {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
42 43 44 45
    // augment time step counter
    ++k;

    // apply finite volume scheme
46
    evolve(grid,mapper,c,t,dt);
Robert Klöfkorn's avatar
Robert Klöfkorn committed
47 48

    // augment time
49
    t += dt;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
50 51 52 53 54

    // check if data should be written
    if (t >= saveStep)
    {
      // write data
55
      vtkout(grid,c,"concentration",counter,t);     /*@\label{fvc:file}@*/
Robert Klöfkorn's avatar
Robert Klöfkorn committed
56 57 58 59 60 61 62

      // increase counter and saveStep for next interval
      saveStep += saveInterval;
      ++counter;
    }

    // print info about time, timestep size and counter
63 64
    std::cout << "s=" << grid.size(0)
              << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
65
  }                                                    /*@\label{fvc:loop1}@*/
66 67 68 69 70 71 72 73
}

//===============================================================
// The main function creates objects and does the time loop
//===============================================================

int main (int argc , char ** argv)
{
Robert Klöfkorn's avatar
Robert Klöfkorn committed
74 75
  // initialize MPI, finalize is done automatically on exit
  Dune::MPIHelper::instance(argc,argv);
76

77 78
  // start try/catch block to get error messages from dune
  try {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
79 80
    using namespace Dune;

81 82 83 84
    // the GridSelector :: GridType is defined in gridtype.hh and is
    // set during compilation
    typedef GridSelector :: GridType Grid;

85
    // use unitcube from dgf grids
Robert Klöfkorn's avatar
Robert Klöfkorn committed
86
    std::stringstream dgfFileName;
87 88
    dgfFileName << DUNE_GRID_HOWTO_EXAMPLE_GRIDS_PATH
      << "unitcube" << Grid::dimension << ".dgf";
Robert Klöfkorn's avatar
Robert Klöfkorn committed
89

90 91
    // create grid pointer
    GridPtr<Grid> gridPtr( dgfFileName.str() );
Robert Klöfkorn's avatar
Robert Klöfkorn committed
92 93

    // grid reference
94
    Grid& grid = *gridPtr;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
95 96

    // half grid width 4 times
97
    int level = 4 * DGFGridInfo<Grid>::refineStepsForHalf();
Robert Klöfkorn's avatar
Robert Klöfkorn committed
98 99 100 101 102 103

    // refine grid until upper limit of level
    grid.globalRefine(level);

    // do time loop until end time 0.5
    timeloop(grid, 0.5);
104 105
  }
  catch (std::exception & e) {
106
    std::cout << "ERROR: " << e.what() << std::endl;
107 108 109 110 111 112
    return 1;
  }
  catch (...) {
    std::cout << "Unknown ERROR" << std::endl;
    return 1;
  }
113 114 115 116

  // done
  return 0;
}