[dune-pdelab] Finite Volume Method with slope reconstruction

[Christian Engwer]

Hi Gregor,

what you are trying to do, is apply a linear operator on your current
solution. You can simply do this by writing a local operator and calling
the residual method on the gridoperator. This one would map from
scalar cell data (potential) to vector valued cell data (slope).

You could alternatively write the assembly yourself, but it would
require exposing some internal implementation details.

Ciao
Christian


On Fri, Apr 28, 2017 at 03:13:02PM +0200, Gregor Corbin wrote:
> Hello dear Dune developers,
> 
> I have a dune-pdelab code that implements a first-order finite volume method
> for a hyperbolic system. Currently I am trying to implement a slope
> reconstruction to achieve second order accuracy in space. Therefore, before
> each stage of the time-stepping scheme, I need to iterate over the grid and
> compute the slopes on each cell.
> My current approach is to write a limiter class and give it to the explicit
> one step method like this:
> 
> oneStepMethod.apply(t,dt,x_old,x_new,LinearReconstruction);
> 
> In the prestage method of this class, which gets only the solution vector, i
> iterate over all the cells on the grid. For each cell I want to extract the
> solution on that cell and on it's neighbors.
> But I have no idea how to do this.
> 
> I appreciate any help on this problem. Is there maybe even a better way to
> implement slope limiting in pdelab?
> 
> Cheers and thanks in advance,
> Gregor Corbin
> 
> Below is a sketch for the code of the limiter class:
> 
> template <typename GridView, typename GFS, int n_components>
> class LinearReconstruction
> {
>   enum {dim = GridView::dimension};
>   enum {n_comp = n_components};
> 
> public:
>   LinearReconstruction(const GridView &grid_view, const GFS &gfs):
>   grid_view_(grid_view),
>   gfs_(gfs),
>   mapper_(grid_view_),
>   slopes_(mapper_.size(),0.0)
>   {}
> 
>   template<typename X>
>   void prestage(X& x)
>   {
>     // reconstruct slopes
>     for (const auto& cell : elements(grid_view_) )
>     {
>        // I would like to do something like in the local operator methods,
>        // but how do I get the lfsu_s ???
>        // is there a method to get it from the cell and the grid function
> space?
> 
>        // lfsu_s = gfs_.localFunctionSpace(cell); <- does something similar
> exist?
>       Dune::FieldVector<double,n_comp> u_s(0.0);
>       for (int k = 0; k < n_comp; ++k)
>         u_s[k] = x(lfsu_s.child(k),0);
> 
>       // get the solutions from all the neighbors
>       std::vector<Dune::FieldVector<double,n_comp> > u_n(pow(2,dim),0.0); //
> only valid for structured cubic grids
>       int i = 0;
>       for (const auto& is : intersections(grid_view_,cell))
>       {
>         const auto &neighbor = is.outside();
>          // lfsu_n = gfs_.localFunctionSpace(neighbor); <- does something
> similar exist?
> 
>          for (int k = 0; k < n_comp; ++k)
>               u_n.at(i)[k] = x(lfsu_n.child(k),0);
>          ++i;
>       }
> 
>       // use the info to reconstruct slopes
>       slopes_.at(mapper_.index(cell)) =  // some function(u_s, u_n)
>     }
>   }
> 
>   template<typename V>
>   void poststage(V& v)
>   {}
> private:
>   const GridView &grid_view_;
>    const GFS &gfs_;
>   using Mapper =
> Dune::MultipleCodimMultipleGeomTypeMapper<GridView,Dune::MCMGElementLayout>;
>   const Mapper mapper_;
>   std::vector<Dune::FieldMatrix<double,n_comp,dim> > slopes_;
> };
> 
> 
> 
> 
> 
> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab
>