parfinitevolume.cc 3.96 KB
Newer Older
Peter Bastian's avatar
Peter Bastian committed
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
7
#include <dune/grid/common/mcmgmapper.hh> // mapper class
Robert Klöfkorn's avatar
Robert Klöfkorn committed
8
#include <dune/common/mpihelper.hh> // include mpi helper class
Peter Bastian's avatar
Peter Bastian committed
9

10 11

// checks for defined gridtype and inlcudes appropriate dgfparser implementation
Robert Klöfkorn's avatar
Robert Klöfkorn committed
12
#include "vtkout.hh"
Peter Bastian's avatar
Peter Bastian committed
13
#include "unitcube.hh"
14
#include "transportproblem2.hh"
Peter Bastian's avatar
Peter Bastian committed
15 16 17 18
#include "initialize.hh"
#include "parfvdatahandle.hh"
#include "parevolve.hh"

Robert Klöfkorn's avatar
Robert Klöfkorn committed
19

Peter Bastian's avatar
Peter Bastian committed
20 21 22 23 24 25 26 27
//===============================================================
// the time loop function working for all types of grids
//===============================================================

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

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

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

  // initialize concentration with initial values
  initialize(grid,mapper,c);
47
  vtkout(grid,c,"pconc",0,0.0,grid.comm().rank());
Peter Bastian's avatar
Peter Bastian committed
48 49 50 51

  // now do the time steps
  double t=0,dt;
  int k=0;
52 53 54
  const double saveInterval = 0.1;
  double saveStep = 0.1;
  int counter = 1;
Peter Bastian's avatar
Peter Bastian committed
55 56
  while (t<tend)
  {
57
    // augment time step counter
Peter Bastian's avatar
Peter Bastian committed
58
    k++;
59 60

    // apply finite volume scheme
Peter Bastian's avatar
Peter Bastian committed
61
    parevolve(grid,mapper,c,t,dt);
62 63

    // augment time
Peter Bastian's avatar
Peter Bastian committed
64
    t += dt;
65 66 67 68 69 70 71 72 73 74 75 76 77

    // check if data should be written
    if (t >= saveStep)
    {
      // write data
      vtkout(grid,c,"pconc",counter,t,grid.comm().rank());

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

    // print info about time, timestep size and counter
78
    if (grid.comm().rank()==0)                         /*@\label{pfc:rank0}@*/
79
      std::cout << "k=" << k << " t=" << t << " dt=" << dt << std::endl;
Peter Bastian's avatar
Peter Bastian committed
80
  }
81
  vtkout(grid,c,"pconc",counter,tend,grid.comm().rank());
Peter Bastian's avatar
Peter Bastian committed
82 83 84 85 86 87 88 89
}

//===============================================================
// 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
90 91
  // initialize MPI, finalize is done automatically on exit
  Dune::MPIHelper::instance(argc,argv);
Peter Bastian's avatar
Peter Bastian committed
92

93 94
  // start try/catch block to get error messages from dune
  try {
95 96 97
    using namespace Dune;

    UnitCube<YaspGrid<2>,64> uc;
98 99
    uc.grid().globalRefine(2);
    partimeloop(uc.grid(),0.5);
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129

    /* To use an alternative grid implementations for parallel computations,
       uncomment exactly one definition of uc2 and the line below. */
    //    #define LOAD_BALANCING

    //  UGGrid supports parallelization in 2 or 3 dimensions
#if HAVE_UG
    //    typedef UGGrid< 2 > GridType;
    //    UnitCube< GridType, 2 > uc2;
#endif

    //  ALUGRID supports parallelization in 3 dimensions only
#if HAVE_ALUGRID
    //    typedef ALUCubeGrid< 3, 3 > GridType;
    //    typedef ALUSimplexGrid< 3, 3 > GridType;
    //    UnitCube< GridType , 1 > uc2;
#endif

#ifdef LOAD_BALANCING

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

    // re-partition grid for better load balancing
    uc2.grid().loadBalance();                               /*@\label{pfv:lb}@*/

    // do time loop until end time 0.5
    partimeloop(uc2.grid(), 0.5);
#endif

130 131 132 133 134 135 136 137 138 139 140 141 142
  }
  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;
  }
Peter Bastian's avatar
Peter Bastian committed
143 144 145 146

  // done
  return 0;
}