[dune-pdelab] Finite Volume Method with slope reconstruction

[Gregor Corbin]

Hi Christian,

thanks for the quick reply.

If I use a grid operator for this, how do I get the solution on 
neighboring cells?  I suppose I could use the alpha_skeleton method to 
communicate this information and then doing the reconstruction in the 
alpha_volume_post_skeleton method. But this would mean a lot of storage 
and copying overhead.

Also, the reconstruction has to be done before each stage of the time 
stepping. So would I have to call the reconstruction operator from the 
preStage() method of my original local operator?

Another question related to this: Are the intersections in a YASPGrid 
ordered in a consistent way? I.e. has intersections(grid_view,cell)[0] 
always the same normal?

Cheers,

Gregor

On 28.04.2017 15:32, Christian Engwer wrote:
> 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
>>