parevolve.hh 5.25 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 4 5 6 7 8 9
#ifndef __DUNE_GRID_HOWTO_PAREVOLVE_HH__
#define __DUNE_GRID_HOWTO_PAREVOLVE_HH__

#include "parfvdatahandle.hh"

#include <dune/grid/common/gridenums.hh>
#include <dune/common/fvector.hh>
10

Peter Bastian's avatar
Peter Bastian committed
11 12 13 14
template<class G, class M, class V>
void parevolve (const G& grid, const M& mapper, V& c, double t, double& dt)
{
  // check data partitioning
15
  assert(grid.overlapSize(0)>0 || (grid.ghostSize(0)>0)); /*@\label{peh:assert}@*/
Peter Bastian's avatar
Peter Bastian committed
16 17 18 19 20 21 22 23

  // first we extract the dimensions of the grid
  const int dim = G::dimension;
  const int dimworld = G::dimensionworld;

  // type used for coordinates in the grid
  typedef typename G::ctype ct;

24 25 26
  // type for grid view on leaf part
  typedef typename G::LeafGridView GridView;

Peter Bastian's avatar
Peter Bastian committed
27
  // iterator type
28 29
  typedef typename GridView::template Codim<0>::
  template Partition<Dune::All_Partition>::Iterator LeafIterator;       /*@\label{peh:pit}@*/
Peter Bastian's avatar
Peter Bastian committed
30

31 32 33
  // leaf entity geometry
  typedef typename LeafIterator::Entity::Geometry LeafGeometry;

Peter Bastian's avatar
Peter Bastian committed
34
  // intersection iterator type
35
  typedef typename GridView::IntersectionIterator IntersectionIterator;
Peter Bastian's avatar
Peter Bastian committed
36

37 38 39
  // type of intersection
  typedef typename IntersectionIterator::Intersection Intersection;

40 41 42
  // intersection geometry
  typedef typename Intersection::Geometry IntersectionGeometry;

43 44
  // entity type
  typedef typename G::template Codim<0>::Entity Entity;
Peter Bastian's avatar
Peter Bastian committed
45 46

  // allocate a temporary vector for the update
Peter Bastian's avatar
Peter Bastian committed
47
  V update(c.size());
Christian Engwer's avatar
Christian Engwer committed
48
  for (typename V::size_type i=0; i<c.size(); i++) update[i] = 0;
Peter Bastian's avatar
Peter Bastian committed
49 50 51 52

  // initialize dt very large
  dt = 1E100;

53
  // get grid view instance on leaf grid
54
  GridView gridView = grid.leafGridView();
55

Peter Bastian's avatar
Peter Bastian committed
56 57
  // compute update vector and optimum dt in one grid traversal
  // iterate over all entities, but update is only used on interior entities
58 59
  LeafIterator endit = gridView.template end<0,Dune::All_Partition>(); /*@\label{peh:end}@*/
  for (LeafIterator it = gridView.template begin<0,Dune::All_Partition>(); it!=endit; ++it) /*@\label{peh:begin}@*/
Peter Bastian's avatar
Peter Bastian committed
60
  {
61
    // cell geometry
62
    const LeafGeometry geo = it->geometry();
Peter Bastian's avatar
Peter Bastian committed
63

64
    // cell volume
65
    double volume = geo.volume();
Peter Bastian's avatar
Peter Bastian committed
66 67

    // cell index
68
    int indexi = mapper.index(*it);
Peter Bastian's avatar
Peter Bastian committed
69 70 71 72 73

    // variable to compute sum of positive factors
    double sumfactor = 0.0;

    // run through all intersections with neighbors and boundary
74 75
    const IntersectionIterator isend = gridView.iend(*it);
    for( IntersectionIterator is = gridView.ibegin(*it); is != isend; ++is )
Peter Bastian's avatar
Peter Bastian committed
76
    {
77 78
      const Intersection &intersection = *is;

Peter Bastian's avatar
Peter Bastian committed
79
      // get geometry type of face
80
      const IntersectionGeometry igeo = intersection.geometry();
Peter Bastian's avatar
Peter Bastian committed
81

82
      // get normal vector scaled with volume
83
      Dune::FieldVector< ct, dimworld > integrationOuterNormal
84
        = intersection.centerUnitOuterNormal();
85
      integrationOuterNormal *= igeo.volume();
Peter Bastian's avatar
Peter Bastian committed
86 87

      // center of face in global coordinates
88
      Dune::FieldVector< ct, dimworld > faceglobal = igeo.center();
Peter Bastian's avatar
Peter Bastian committed
89 90 91 92 93 94 95 96 97 98 99

      // evaluate velocity at face center
      Dune::FieldVector<double,dim> velocity = u(faceglobal,t);

      // compute factor occuring in flux formula
      double factor = velocity*integrationOuterNormal/volume;

      // for time step calculation
      if (factor>=0) sumfactor += factor;

      // handle interior face
100
      if( intersection.neighbor() )
Peter Bastian's avatar
Peter Bastian committed
101 102
      {
        // access neighbor
103 104
        Entity outside = intersection.outside();
        int indexj = mapper.index(outside);
Peter Bastian's avatar
Peter Bastian committed
105

106
        const int insideLevel = it->level();
107
        const int outsideLevel = outside.level();
108

109
        // handle face from one side
110 111
        if( (insideLevel > outsideLevel)
            || ((insideLevel == outsideLevel) && (indexi < indexj)) )
Peter Bastian's avatar
Peter Bastian committed
112
        {
113
          // compute factor in neighbor
114
          const LeafGeometry nbgeo = outside.geometry();
115
          double nbvolume = nbgeo.volume();
116 117
          double nbfactor = velocity*integrationOuterNormal/nbvolume;

118
          if( factor < 0 )       // inflow
Peter Bastian's avatar
Peter Bastian committed
119
          {
120 121 122
            update[indexi] -= c[indexj]*factor;
            update[indexj] += c[indexj]*nbfactor;
          }
123
          else       // outflow
124 125 126
          {
            update[indexi] -= c[indexi]*factor;
            update[indexj] += c[indexi]*nbfactor;
Peter Bastian's avatar
Peter Bastian committed
127 128 129 130 131
          }
        }
      }

      // handle boundary face
132 133 134
      if( intersection.boundary() )
      {
        if( factor < 0 )       // inflow, apply boundary condition
Peter Bastian's avatar
Peter Bastian committed
135
          update[indexi] -= b(faceglobal,t)*factor;
136
        else       // outflow
Peter Bastian's avatar
Peter Bastian committed
137
          update[indexi] -= c[indexi]*factor;
138 139
      }
    }       // end all intersections
Peter Bastian's avatar
Peter Bastian committed
140 141

    // compute dt restriction
142
    if (it->partitionType()==Dune::InteriorEntity)       /*@\label{peh:inter}@*/
Peter Bastian's avatar
Peter Bastian committed
143 144 145 146
      dt = std::min(dt,1.0/sumfactor);

  }       // end grid traversal

147 148
  // global min over all partitions
  dt = grid.comm().min(dt);                            /*@\label{peh:min}@*/
Peter Bastian's avatar
Peter Bastian committed
149 150 151 152
  // scale dt with safety factor
  dt *= 0.99;

  // exchange update
153
  VectorExchange<M,V> dh(mapper,update);               /*@\label{peh:dist0}@*/
154 155
  grid.template
  communicate<VectorExchange<M,V> >(dh,Dune::InteriorBorder_All_Interface,
156
                                    Dune::ForwardCommunication);     /*@\label{peh:dist1}@*/
Peter Bastian's avatar
Peter Bastian committed
157 158 159 160 161 162 163

  // update the concentration vector
  for (unsigned int i=0; i<c.size(); ++i)
    c[i] += dt*update[i];

  return;
}
164 165

#endif //__DUNE_GRID_HOWTO_PAREVOLVE_HH__