fokkerplanck.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#include <memory>
#include <vector>
#include <array>
#include <random>
#include <cmath>
#include <iostream>
 
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
 
 
#include <dumux/parallel/parallel_for.hh>
 
#include <dumux/geometry/intersectingentities.hh>
#include <dumux/geometry/geometryintersection.hh>
 
#include <dumux/particles/particle.hh>
#include <dumux/particles/simpleparticlecloud.hh>
 
namespace Dumux::Particles {
 
template<class Scalar, class GridGeometry>
class FokkerPlanck
{
    static constexpr int dim = GridGeometry::GridView::dimension;
    static constexpr int dimWorld = GridGeometry::GridView::dimensionworld;
 
    using Particle = Dumux::Particles::Particle<dimWorld, Scalar>;
    using Cloud = SimpleParticleCloud<Particle>;
    using GlobalPosition = Particle::GlobalPosition;
    using Velocity = Dune::FieldVector<Scalar, dim>;
 
    static_assert(GridGeometry::discMethod == DiscretizationMethods::cctpfa, "Implementation currently assumes a CC-TPFA grid geometry");
public:
    FokkerPlanck(std::shared_ptr<const GridGeometry> gg, std::shared_ptr<Cloud> cloud)
    : gridGeometry_(std::move(gg))
    , cloud_(std::move(cloud))
    , rndGen_(std::random_device{}())
    {
        if (cloud_->size(false) == 0)
            cloud_->resize(getParam<int>("FokkerPlanck.NumParticles"));
 
        const auto numParticles = cloud_->size(false);
        diffusionCoefficient_ = getParam<Scalar>("FokkerPlanck.DiffusionCoefficient");
        data_.resize(numParticles);
        particleMass_ = getParam<Scalar>("FokkerPlanck.ParticleMass", 1.0/numParticles);
        volumes_.resize(gridGeometry_->gridView().size(0), 0.0);
        for (const auto& element : elements(gridGeometry_->gridView()))
        {
            const auto eIdx = gridGeometry_->elementMapper().index(element);
            const auto geometry = element.geometry();
            volumes_[eIdx] = geometry.volume();
        }
    }
 
    struct ParticleData
    {
        std::size_t eIdx; // in which element of the grid this particle is
        Scalar time; // at which time the particle lives
    };
 
    template<class Func>
    void init(const Func& func)
    {
        cloud_->forEachParticle([&](auto& p){
            func(p, data_[p.id()]);
        });
    }
 
    void run(const Scalar dt)
    {
        if (velocities_.empty())
            DUNE_THROW(Dune::InvalidStateException, "A velocity field must be set before calling run");
 
        // move the particles with the given velocity field through the grid
        cloud_->forEachActiveParticle([&](auto& p){
            moveParticle_(p, dt);
        });
    }
 
    void computeDensity(std::vector<Scalar>& density) const
    {
        density.assign(gridGeometry_->gridView().size(0), 0.0);
        for (const auto& p : particles(*cloud_))
        {
            const auto eIdx = data_[p.id()].eIdx;
            density[eIdx] += particleMass_/volumes_[eIdx];
        }
    }
 
    void setVelocityData(const std::vector<Velocity>& velocities)
    {
        if (velocities.size() != gridGeometry_->gridView().size(0) && velocities.size() != 1)
            DUNE_THROW(Dune::InvalidStateException, "Velocity data size does not match number of elements or 1");
        velocities_ = velocities;
    }
 
    std::shared_ptr<const Cloud> particleCloud() const
    { return cloud_; }
 
private:
    void moveParticle_(Particle& p, const Scalar dt)
    {
        auto& particleData = data_[p.id()];
 
        // (1) Drift: time loop until time integration is finished
        particleData.time = dt; // recharge time
        int count = 0;
        while (particleData.time > 0.0 && velocity_(particleData.eIdx).two_norm2() > 0.0)
        {
            const auto tStartIteration = particleData.time;
 
            const auto velocityInGrid = velocity_(particleData.eIdx);
            moveParticleInGridElement_(p, velocityInGrid);
 
            count = particleData.time < tStartIteration ? 0 : count + 1;
            if (count > 4)
            {
                std::cout << "Particle [" << p.id() << "] at " << p.position() << " in element " << particleData.eIdx << ", t=" << particleData.time << ", velocity: " << velocityInGrid << "\n";
                std::cout << "-- during advection step: particle got stuck and is deactivated\n";
 
                cloud_->deactivateParticle(p);
                particleData.time = -1.0;
            }
        }
 
        static_assert(dim == dimWorld, "Not implemented yet for dim != dimWorld");
 
        // (2) Brownian motion: time loop until time integration is finished
        particleData.time = 1.0; // recharge time
        const auto velocityInGrid = [&]{
            const auto scale = std::sqrt(2.0*diffusionCoefficient_);
            Dune::FieldVector<Scalar, dimWorld> newPos(p.position());
            for (int i = 0; i < dim; ++i)
                newPos[i] += scale * std::normal_distribution<Scalar>(0.0, std::sqrt(dt))(rndGen_);
            return newPos - p.position();
        }();
 
        const auto velocityNorm = velocityInGrid.two_norm2();
        count = 0;
        while (particleData.time > 0.0 && velocityNorm > 0.0)
        {
            const auto tStartIteration = particleData.time;
 
            moveParticleInGridElement_(p, velocityInGrid);
 
            count = particleData.time < tStartIteration ? 0 : count + 1;
            if (count > 4)
            {
                std::cout << "Particle [" << p.id() << "] at " << p.position() << " in element " << particleData.eIdx << ", t=" << particleData.time << ", velocity: " << velocityInGrid << "\n";
                std::cout << "-- during diffusion step: particle got stuck and is deactivated\n";
 
                cloud_->deactivateParticle(p);
                particleData.time = -1.0;
            }
        }
    }
 
    void moveParticleInGridElement_(Particle& p, const Velocity& velocity)
    {
        auto& particleData = data_[p.id()];
        assert(particleData.time > 0.0);
 
        const auto eIdx = particleData.eIdx;
        if (eIdx >= gridGeometry_->gridView().size(0))
            DUNE_THROW(Dune::InvalidStateException, "Element index is invalid");
 
        const auto& element = gridGeometry_->element(eIdx);
        const auto geometry = element.geometry();
 
        const auto startPos = p.position();
        auto endPos = startPos;
        endPos.axpy(particleData.time, velocity);
 
        // we stay inside this element
        if (intersectsPointGeometry(endPos, geometry))
        {
            p.setPosition(endPos);
            particleData.time = -1.0; // this particle is done
            return;
        }
 
        // increase the change of finding the intersection
        // if we are right on the element facet by backtracking a bit
        auto beforeStart = startPos;
        beforeStart.axpy(-1e-4*particleData.time, velocity);
 
        const auto fvGeometry = localView(*gridGeometry_).bind(element);
        Dune::AffineGeometry<double, 1, dimWorld> segGeo(Dune::GeometryTypes::line, std::array<GlobalPosition, 2>{{beforeStart, endPos}});
        using IP = Dumux::IntersectionPolicy::PointPolicy<double, dimWorld>;
        using IntersectionAlgorithm = GeometryIntersection<Dune::AffineGeometry<double, 1, dimWorld>, decltype(fvGeometry.geometry(fvGeometry.scvf(0))), IP>;
        using Intersection = typename IntersectionAlgorithm::Intersection;
        Intersection intersection;
 
        static_assert(dim >= 1, "Particle transport only makes sense for 1D and higher");
        const auto intersectionFound = [&](const auto& scvf, Intersection& intersection)
        {
            if constexpr (dim == 1)
            {
                const bool found = intersectsPointGeometry(scvf.ipGlobal(), segGeo);
                if (found) intersection = scvf.ipGlobal();
                return found;
            }
            else
                return IntersectionAlgorithm::intersection(
                    segGeo, fvGeometry.geometry(scvf), intersection
                );
        };
 
        for (const auto& scvf : scvfs(fvGeometry))
        {
            if (intersectionFound(scvf, intersection))
            {
                // for now we deactivate the particle if it hits a boundary
                if (scvf.boundary())
                {
                    cloud_->deactivateParticle(p);
                    particleData.time = -1.0; // this particle is done
                    return;
                }
 
                // corner case (intersects more than one scvfs)
                if (scvf.unitOuterNormal()*velocity < 1e-4*velocity.two_norm())
                    continue;
 
                // put particle on face into the neighboring element
                const auto velocityNorm = velocity.two_norm2();
                const auto distance = (intersection-startPos).two_norm2();
                intersection.axpy(1e-3, endPos-startPos);
                p.setPosition(intersection);
                particleData.eIdx = scvf.outsideScvIdx();
                particleData.time -= std::sqrt(distance/velocityNorm);
                return;
            }
        }
 
        // particle got lost: try to find the particle in the grid
        static constexpr bool isCartesianGrid = Dune::Capabilities::isCartesian<typename GridGeometry::GridView::Grid>::v;
        const auto entities = intersectingEntities(p.position(), gridGeometry_->boundingBoxTree(), isCartesianGrid);
        if (entities.empty()) // outside the grid
        {
            cloud_->deactivateParticle(p);
            particleData.time = -1.0; // this particle is done
            return;
        }
 
        // found particle in an element, try to move on from there
        particleData.eIdx = entities[0];
    }
 
    const Velocity& velocity_(std::size_t eIdx) const
    {
        if (velocities_.size() == 1)
            return velocities_[0];
        return velocities_[eIdx];
    }
 
    std::shared_ptr<const GridGeometry> gridGeometry_;
    std::vector<ParticleData> data_;
    std::shared_ptr<Cloud> cloud_;
    std::mt19937 rndGen_;
 
    std::vector<Velocity> velocities_;
    Scalar diffusionCoefficient_;
 
    std::vector<Scalar> volumes_; 
    Scalar particleMass_;
 
};
 
} // end namespace Dumux::Particles