parfinitevolume.cc 3.78 KB
Newer Older
Peter Bastian's avatar
Peter Bastian committed
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
Peter Bastian's avatar
Peter Bastian committed
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
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
//===============================================================

template<class G>
void partimeloop (const G& grid, double tend)
{
  // make a mapper for codim 0 entities in the leaf grid
28 29
  Dune::LeafMultipleCodimMultipleGeomTypeMapper<G>
  mapper(grid, Dune::mcmgElementLayout());
Peter Bastian's avatar
Peter Bastian committed
30 31 32 33 34 35

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

  // initialize concentration with initial values
  initialize(grid,mapper,c);
36
  vtkout(grid,c,"pconc",0,0.0,grid.comm().rank());
Peter Bastian's avatar
Peter Bastian committed
37 38 39 40

  // now do the time steps
  double t=0,dt;
  int k=0;
41 42 43
  const double saveInterval = 0.1;
  double saveStep = 0.1;
  int counter = 1;
Peter Bastian's avatar
Peter Bastian committed
44 45
  while (t<tend)
  {
46
    // augment time step counter
Peter Bastian's avatar
Peter Bastian committed
47
    k++;
48 49

    // apply finite volume scheme
Peter Bastian's avatar
Peter Bastian committed
50
    parevolve(grid,mapper,c,t,dt);
51 52

    // augment time
Peter Bastian's avatar
Peter Bastian committed
53
    t += dt;
54 55 56 57 58 59 60 61 62 63 64 65 66

    // 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
67
    if (grid.comm().rank()==0)                         /*@\label{pfc:rank0}@*/
68
      std::cout << "k=" << k << " t=" << t << " dt=" << dt << std::endl;
Peter Bastian's avatar
Peter Bastian committed
69
  }
70
  vtkout(grid,c,"pconc",counter,tend,grid.comm().rank());
Peter Bastian's avatar
Peter Bastian committed
71 72 73 74 75 76 77 78
}

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

82 83
  // start try/catch block to get error messages from dune
  try {
84 85 86
    using namespace Dune;

    UnitCube<YaspGrid<2>,64> uc;
87 88
    uc.grid().globalRefine(2);
    partimeloop(uc.grid(),0.5);
89 90 91 92 93 94 95 96 97 98 99

    /* 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

100
    //  ALUGRID supports parallelization in 2 or 3 dimensions
101
#if HAVE_DUNE_ALUGRID
102 103
    //    typedef Dune::ALUGrid< 3, 3, Dune::cube, Dune::nonconforming > GridType;
    //    typedef Dune::ALUGrid< 3, 3, Dune::simplex, Dune::nonconforming > GridType;
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
    //    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

119 120
  }
  catch (std::exception & e) {
121
    std::cout << "ERROR: " << e.what() << std::endl;
122 123 124 125 126 127
    return 1;
  }
  catch (...) {
    std::cout << "Unknown ERROR" << std::endl;
    return 1;
  }
Peter Bastian's avatar
Peter Bastian committed
128 129 130 131

  // done
  return 0;
}