finitevolume.cc 3.36 KB
 Peter Bastian committed Mar 17, 2006 1 2 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- // vi: set et ts=4 sw=2 sts=2:  Christoph Grüninger committed Oct 24, 2014 3 #include // know what grids are present  Peter Bastian committed Mar 17, 2006 4 5 6 #include // for input/output to shell #include // for input/output to files #include // STL vector class  Christian Engwer committed May 16, 2006 7 #include // mapper class  Christian Engwer committed Nov 06, 2012 8 #include // include mpi helper class  Robert Klöfkorn committed Jul 29, 2007 9 10  #include "vtkout.hh"  Peter Bastian committed Mar 28, 2006 11 #include "transportproblem2.hh"  Peter Bastian committed Mar 17, 2006 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 void timeloop (const G& grid, double tend) { // make a mapper for codim 0 entities in the leaf grid  Ansgar Burchardt committed Aug 09, 2017 23 24  Dune::LeafMultipleCodimMultipleGeomTypeMapper mapper(grid, Dune::mcmgElementLayout());  Peter Bastian committed Mar 17, 2006 25 26 27 28 29  // allocate a vector for the concentration std::vector c(mapper.size()); // initialize concentration with initial values  Sven Marnach committed Nov 10, 2006 30  initialize(grid,mapper,c); /*@\label{fvc:init}@*/  Christian Engwer committed Mar 09, 2009 31  vtkout(grid,c,"concentration",0,0.0);  Peter Bastian committed Mar 17, 2006 32 33 34 35  // now do the time steps double t=0,dt; int k=0;  Robert Klöfkorn committed Jul 29, 2007 36 37  const double saveInterval = 0.1; double saveStep = 0.1;  Christian Engwer committed Mar 09, 2009 38  int counter = 1;  Robert Klöfkorn committed Jul 29, 2007 39   Sven Marnach committed Nov 10, 2006 40  while (t= saveStep) { // write data  Oliver Sander committed Aug 19, 2014 55  vtkout(grid,c,"concentration",counter,t); /*@\label{fvc:file}@*/  Robert Klöfkorn committed Jul 29, 2007 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  Christian Engwer committed Mar 09, 2009 63 64  std::cout << "s=" << grid.size(0) << " k=" << k << " t=" << t << " dt=" << dt << std::endl;  Sven Marnach committed Nov 10, 2006 65  } /*@\label{fvc:loop1}@*/  Peter Bastian committed Mar 17, 2006 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 committed Jul 29, 2007 74 75  // initialize MPI, finalize is done automatically on exit Dune::MPIHelper::instance(argc,argv);  Peter Bastian committed Mar 17, 2006 76   Christian Engwer committed Dec 05, 2006 77 78  // start try/catch block to get error messages from dune try {  Robert Klöfkorn committed Jul 29, 2007 79 80  using namespace Dune;  Andreas Dedner committed Jan 14, 2010 81 82 83 84  // the GridSelector :: GridType is defined in gridtype.hh and is // set during compilation typedef GridSelector :: GridType Grid;  Christian Engwer committed Mar 09, 2009 85  // use unitcube from dgf grids  Robert Klöfkorn committed Jul 29, 2007 86  std::stringstream dgfFileName;  Christoph Grüninger committed May 10, 2013 87 88  dgfFileName << DUNE_GRID_HOWTO_EXAMPLE_GRIDS_PATH << "unitcube" << Grid::dimension << ".dgf";  Robert Klöfkorn committed Jul 29, 2007 89   Andreas Dedner committed Jan 14, 2010 90 91  // create grid pointer GridPtr gridPtr( dgfFileName.str() );  Robert Klöfkorn committed Jul 29, 2007 92 93  // grid reference  Andreas Dedner committed Jan 14, 2010 94  Grid& grid = *gridPtr;  Robert Klöfkorn committed Jul 29, 2007 95 96  // half grid width 4 times  Andreas Dedner committed Jan 14, 2010 97  int level = 4 * DGFGridInfo::refineStepsForHalf();  Robert Klöfkorn committed Jul 29, 2007 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);  Christian Engwer committed Dec 05, 2006 104 105  } catch (std::exception & e) {  Christoph Grüninger committed Mar 08, 2016 106  std::cout << "ERROR: " << e.what() << std::endl;  Christian Engwer committed Dec 05, 2006 107 108 109 110 111 112  return 1; } catch (...) { std::cout << "Unknown ERROR" << std::endl; return 1; }  Peter Bastian committed Mar 17, 2006 113 114 115 116  // done return 0; }