adaptivefinitevolume.cc 4.32 KB
Newer Older
1 2 3 4 5 6
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"               // know what grids are present
#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

// checks for defined gridtype and inlcudes appropriate dgfparser implementation
Robert Klöfkorn's avatar
Robert Klöfkorn committed
9
#include <dune/grid/io/file/dgfparser/dgfgridtype.hh>
Robert Klöfkorn's avatar
Robert Klöfkorn committed
10

11
#include <dune/grid/common/mcmgmapper.hh> // mapper class
Robert Klöfkorn's avatar
Robert Klöfkorn committed
12
#include <dune/common/mpihelper.hh> // include mpi helper class
13

14
#include "vtkout.hh"
15
#include "unitcube.hh"
16
#include "transportproblem2.hh"
17 18 19 20 21 22 23 24 25 26 27 28
#include "initialize.hh"
#include "evolve.hh"
#include "finitevolumeadapt.hh"

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

//! Parameter for mapper class
template<int dim>
struct P0Layout
{
29
  bool contains (Dune::GeometryType gt)
30
  {
31
    if (gt.dim()==dim) return true;
32 33 34 35 36 37 38 39
    return false;
  }
};

template<class G>
void timeloop (G& grid, double tend, int lmin, int lmax)
{
  // make a mapper for codim 0 entities in the leaf grid
Peter Bastian's avatar
Peter Bastian committed
40 41
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G,P0Layout>
  mapper(grid);
42 43 44 45 46 47 48 49

  // 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
50
    if (grid.maxLevel()>=lmax) break;
51
    finitevolumeadapt(grid,mapper,c,lmin,lmax,0);           /*@\label{afv:in}@*/
52 53
    initialize(grid,mapper,c);
  }
Robert Klöfkorn's avatar
Robert Klöfkorn committed
54 55

  // write initial data
56
  vtkout(grid,c,"concentration",0,0);
57

Robert Klöfkorn's avatar
Robert Klöfkorn committed
58
  // variables for time, timestep etc.
59
  double dt, t=0;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
60 61 62 63 64 65
  double saveStep = 0.1;
  const double saveInterval = 0.1;
  int counter = 0;
  int k = 0;

  std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << std::endl;
66 67
  while (t<tend)
  {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
68 69 70 71
    // augment time step counter
    ++k;

    // apply finite volume scheme
72
    evolve(grid,mapper,c,t,dt);
Robert Klöfkorn's avatar
Robert Klöfkorn committed
73 74

    // augment time
75
    t += dt;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
76 77 78 79 80

    // check if data should be written
    if (t >= saveStep)
    {
      // write data
81
      vtkout(grid,c,"concentration",counter,t);
Robert Klöfkorn's avatar
Robert Klöfkorn committed
82 83 84 85 86 87 88

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

    // print info about time, timestep size and counter
89 90
    std::cout << "s=" << grid.size(0)
              << " k=" << k << " t=" << t << " dt=" << dt << std::endl;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
91 92 93 94

    // for unstructured grids call adaptation algorithm
    if( Dune :: Capabilities :: IsUnstructured<G> :: v )
    {
95
      finitevolumeadapt(grid,mapper,c,lmin,lmax,k);          /*@\label{afv:ad}@*/
Robert Klöfkorn's avatar
Robert Klöfkorn committed
96
    }
97
  }
Robert Klöfkorn's avatar
Robert Klöfkorn committed
98 99

  // write last time step
100 101 102
  vtkout(grid,c,"concentration",counter,tend);

  // write
103 104 105 106 107 108 109 110
}

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

114 115
  // start try/catch block to get error messages from dune
  try {
Robert Klöfkorn's avatar
Robert Klöfkorn committed
116 117
    using namespace Dune;

118 119 120 121
    // 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
122 123
    // use unitcube from grids
    std::stringstream dgfFileName;
124
    dgfFileName << "grids/unitcube" << Grid :: dimension << ".dgf";
Robert Klöfkorn's avatar
Robert Klöfkorn committed
125

126 127
    // create grid pointer
    GridPtr<Grid> gridPtr( dgfFileName.str() );
Robert Klöfkorn's avatar
Robert Klöfkorn committed
128 129

    // grid reference
130
    Grid& grid = *gridPtr;
Robert Klöfkorn's avatar
Robert Klöfkorn committed
131 132

    // minimal allowed level during refinement
133
    int minLevel = 2 * DGFGridInfo<Grid>::refineStepsForHalf();
Robert Klöfkorn's avatar
Robert Klöfkorn committed
134 135 136 137 138

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

    // maximal allowed level during refinement
139
    int maxLevel = minLevel + 3 * DGFGridInfo<Grid>::refineStepsForHalf();
Robert Klöfkorn's avatar
Robert Klöfkorn committed
140 141 142

    // do time loop until end time 0.5
    timeloop(grid, 0.5, minLevel, maxLevel);
143 144 145 146 147 148 149 150 151 152 153 154 155
  }
  catch (std::exception & e) {
    std::cout << "STL ERROR: " << e.what() << std::endl;
    return 1;
  }
  catch (Dune::Exception & e) {
    std::cout << "DUNE ERROR: " << e.what() << std::endl;
    return 1;
  }
  catch (...) {
    std::cout << "Unknown ERROR" << std::endl;
    return 1;
  }
156 157 158 159

  // done
  return 0;
}