adaptivefinitevolume.cc 3.88 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
Robert Klöfkorn's avatar
Robert Klöfkorn committed
7

8
#include <dune/grid/common/mcmgmapper.hh> // mapper class
9
#include <dune/common/parallel/mpihelper.hh> // include mpi helper class
10

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

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

template<class G>
void timeloop (G& grid, double tend, int lmin, int lmax)
{
  // make a mapper for codim 0 entities in the leaf grid
25 26
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G>
  mapper(grid, Dune::mcmgElementLayout());
27 28 29 30 31 32 33 34

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

  // initialize concentration with initial values
  initialize(grid,mapper,c);
  for (int i=grid.maxLevel(); i<lmax; i++)
  {
Peter Bastian's avatar
Peter Bastian committed
35
    if (grid.maxLevel()>=lmax) break;
36
    finitevolumeadapt(grid,mapper,c,lmin,lmax,0);           /*@\label{afv:in}@*/
37 38
    initialize(grid,mapper,c);
  }
Robert Klöfkorn's avatar
Robert Klöfkorn committed
39 40

  // write initial data
41
  vtkout(grid,c,"concentration",0,0);
42

Robert Klöfkorn's avatar
Robert Klöfkorn committed
43
  // variables for time, timestep etc.
44
  double dt, t=0;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
45
  double saveStep = 0.1;
Andreas Dedner's avatar
Andreas Dedner committed
46 47
  const double saveInterval = t + 0.1;
  int counter = 1;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
48 49 50
  int k = 0;

  std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << std::endl;
51 52
  while (t<tend)
  {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
53 54 55 56
    // augment time step counter
    ++k;

    // apply finite volume scheme
57
    evolve(grid,mapper,c,t,dt);
Robert Klöfkorn's avatar
Robert Klöfkorn committed
58 59

    // augment time
60
    t += dt;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
61 62 63 64 65

    // check if data should be written
    if (t >= saveStep)
    {
      // write data
66
      vtkout(grid,c,"concentration",counter,t);
Robert Klöfkorn's avatar
Robert Klöfkorn committed
67 68 69 70 71 72 73

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

    // print info about time, timestep size and counter
74 75
    std::cout << "s=" << grid.size(0)
              << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
76 77

    // for unstructured grids call adaptation algorithm
78
    finitevolumeadapt(grid,mapper,c,lmin,lmax,k);        /*@\label{afv:ad}@*/
79
  }
Robert Klöfkorn's avatar
Robert Klöfkorn committed
80 81

  // write last time step
82 83 84
  vtkout(grid,c,"concentration",counter,tend);

  // write
85 86 87 88 89 90 91 92
}

//===============================================================
// 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
93 94
  // initialize MPI, finalize is done automatically on exit
  Dune::MPIHelper::instance(argc,argv);
95

96 97
  // start try/catch block to get error messages from dune
  try {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
98 99
    using namespace Dune;

100 101 102 103
    // the GridSelector :: GridType is defined in gridtype.hh and is
    // set during compilation
    typedef GridSelector :: GridType Grid;

Robert Klöfkorn's avatar
Robert Klöfkorn committed
104 105
    // use unitcube from grids
    std::stringstream dgfFileName;
106 107
    dgfFileName << DUNE_GRID_HOWTO_EXAMPLE_GRIDS_PATH
      << "unitcube" << Grid::dimension << ".dgf";
Robert Klöfkorn's avatar
Robert Klöfkorn committed
108

109 110
    // create grid pointer
    GridPtr<Grid> gridPtr( dgfFileName.str() );
Robert Klöfkorn's avatar
Robert Klöfkorn committed
111 112

    // grid reference
113
    Grid& grid = *gridPtr;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
114 115

    // minimal allowed level during refinement
116
    int minLevel = 2 * DGFGridInfo<Grid>::refineStepsForHalf();
Robert Klöfkorn's avatar
Robert Klöfkorn committed
117 118 119 120 121

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

    // maximal allowed level during refinement
122
    int maxLevel = minLevel + 3 * DGFGridInfo<Grid>::refineStepsForHalf();
Robert Klöfkorn's avatar
Robert Klöfkorn committed
123 124 125

    // do time loop until end time 0.5
    timeloop(grid, 0.5, minLevel, maxLevel);
126 127
  }
  catch (std::exception & e) {
128
    std::cout << "ERROR: " << e.what() << std::endl;
129 130 131 132 133 134
    return 1;
  }
  catch (...) {
    std::cout << "Unknown ERROR" << std::endl;
    return 1;
  }
135 136 137 138

  // done
  return 0;
}