FluidSolver.cpp

/* Copyright (c) 2013       David Sichau <mail"at"sichau"dot"eu>
 *               2013-2015  Simon Tanaka <tanakas"at"gmx"dot"ch>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */
#include <LbmLib/include/Constants.hpp>
#include <LbmLib/include/GlobalSimulationParameters.hpp>
#include <LbmLib/include/nodes/PhysicalNode.hpp>
#include <LbmLib/include/solver/FluidSolver/FluidSolver.hpp>
#include <LbmLib/include/solver/FluidSolver/FluidSolver.hpp>
#include <UtilLib/include/Log.hpp>
#include <algorithm>
#include <cassert>
#include <numeric>
#include <sstream>
#include <string>
namespace LbmLib {
namespace solver {
namespace {
constexpr double W0 = 4.0 / 9.0;
constexpr double W1 = 1.0 / 9.0;
constexpr double W2 = 1.0 / 36.0;
}

void FluidSolver::initSolver() {
    if (velocityInit_) {
        initWithVelocity();
    } else {
        // init the populations with rho = 1 and u=0 and v=0
        // hack to set rho to 1 for the init step

        distributions_[T] = 1.0;
        collide();
    }
}

void FluidSolver::initWithVelocity() {
    const double rho = 1.0;

    // calculate the speeds
    const double u = velocity_.x;
    const double v = velocity_.y;

    const double u2 = u * u;
    const double v2 = v * v;

    const double W1rho = W1 * rho;
    const double W2rho = W2 * rho;
    const double W3 = 9.0 / 2.0;
    const double lastTerm = 3.0 / 2.0 * (u2 + v2);

    double temp[9];

    temp[T] = W0 * rho * (1.0 - lastTerm);
    temp[E] = W1rho * (1.0 + 3.0 * u + W3 * u2 - lastTerm);
    temp[NE] = W2rho * (1.0 + 3.0 * (u + v) + W3 * (u + v) * (u + v) - lastTerm);
    temp[N] = W1rho * (1.0 + 3.0 * v + W3 * v2 - lastTerm);
    temp[NW] = W2rho *
        (1.0 + 3.0 * (-u + v) + W3 * (-u + v) * (-u + v) - lastTerm);
    temp[W] = W1rho * (1.0 + 3.0 * (-u) + W3 * u2 - lastTerm);
    temp[SW] = W2rho *
        (1.0 + 3.0 * (-u - v) + W3 * (-u - v) * (-u - v) - lastTerm);
    temp[S] = W1rho * (1.0 + 3.0 * (-v) + W3 * v2 - lastTerm);
    temp[SE] = W2rho * (1.0 + 3.0 * (u - v) + W3 * (u - v) * (u - v) - lastTerm);

    for (auto d : dirIter_) {
        // assign the equilibirum to the distribution
        // make relaxation
        distributions_[d] = temp[d];
    }


    // reseting the rho to add to zero as it is shifted to the distributions
    rhoToAdd_ = 0.0;
    localSwap();
}

void FluidSolver::collide() {

    // Calculate the rho
    double rho = this->getRho();
    assert(rho>0.0);
    if (std::isfinite(this->rhoToAdd_)) {
        rho += this->rhoToAdd_;
    }
    if (rho>0.0) {

    }
    else {
        std::cout<<"rhoToAdd="<<this->rhoToAdd_<<std::endl;
    }
    assert(rho>0.0);
    double rhoI = 1.0 / rho;
    assert(std::isfinite(rhoI));

    // calculate the velocity
    Field<double> u;
    computeVelocity(rho,this->distributions_,this->force_,u);
    const double ux2 = u.x * u.x;
    const double uy2 = u.y * u.y;


    const double lastTerm = 3.0 / 2.0 * (ux2 + uy2);
    const double tauI = 1.0 / getTau();

    double temp[9];

    const double W3 = 9.0 / 2.0;
    const double W1rho = W1 * rho;
    const double W2rho = W2 * rho;

    // the force calculation is added to the temp calculation
    // an if condition to check if the force is zero is the same speed as simply calculate this
    // A.A. Mohamad, A. Kuzmin (2010): A critical evaluation of force term in lattice Boltzmann method, natural convection problem
    std::vector<double> forcecontribution;
    this->computeForce(this->getTau(),
                       this->getRho(),
                       this->getVelocity(),
                       this->force_,
                       forcecontribution);

    temp[T] = W0 * rho * (1.0 - lastTerm) * tauI +
            forcecontribution[T];
    temp[E] = W1rho * (1.0 + 3.0 * u.x + W3 * ux2 - lastTerm) * tauI +
            forcecontribution[E];
    temp[NE] = W2rho * (1.0 + 3.0 * (u.x + u.y) + W3 * (u.x + u.y) * (u.x + u.y) - lastTerm) * tauI +
            forcecontribution[NE];
    temp[N] = W1rho * (1.0 + 3.0 * u.y + W3 * uy2 - lastTerm) * tauI +
            forcecontribution[N];
    temp[NW] = W2rho * (1.0 + 3.0 * (-u.x + u.y) + W3 * (-u.x + u.y) * (-u.x + u.y) - lastTerm) * tauI +
            forcecontribution[NW];
    temp[W] = W1rho * (1.0 + 3.0 * (-u.x) + W3 * ux2 - lastTerm) * tauI +
            forcecontribution[W];
    temp[SW] = W2rho * (1.0 + 3.0 * (-u.x - u.y) + W3 * (-u.x - u.y) * (-u.x - u.y) - lastTerm) * tauI +
            forcecontribution[SW];
    temp[S] = W1rho * (1.0 + 3.0 * (-u.y) + W3 * uy2 - lastTerm) * tauI +
            forcecontribution[S];
    temp[SE] = W2rho * (1.0 + 3.0 * (u.x - u.y) + W3 * (u.x - u.y) * (u.x - u.y) - lastTerm) * tauI +
            forcecontribution[SE];

    for (auto d : dirIter_) {
        const double tempD = distributions_[d];
        // compute non equilibirum and make relaxation
        distributions_[d] = tempD - tempD * tauI + temp[d];
    }

    // recalculate the velocity
    computeVelocity(rho,this->distributions_,this->force_,u);
    assert(std::isfinite(u.x));
    assert(std::isfinite(u.y));
    velocity_.x = u.x;
    velocity_.y = u.y;

    // reseting the rho to add to zero as it is shifted to the distributions
    rhoToAdd_ = 0.0;
    // preparation for advect step
    localSwap();
}

void FluidSolver::advect() {
    std::swap(distributions_[E], physicalNode_.getPhysicalNeighbour(
                    W)->getFluidSolver().accessDistribution(W));
    std::swap(distributions_[NE], physicalNode_.getPhysicalNeighbour(
                    SW)->getFluidSolver().accessDistribution(SW));
    std::swap(distributions_[N], physicalNode_.getPhysicalNeighbour(
                    S)->getFluidSolver().accessDistribution(S));
    std::swap(distributions_[NW], physicalNode_.getPhysicalNeighbour(
                    SE)->getFluidSolver().accessDistribution(SE));
}

void FluidSolver::addForce(Field<double> f) {
    force_ += f;
}

void FluidSolver::writeSolver(std::ostream* const stream) {
    (*stream) << physicalNode_.getXPos() << '\t' << physicalNode_.getYPos();
    for (auto d : distributions_) {
        (*stream) << '\t' << d;
    }
    (*stream) << '\n';
}

void FluidSolver::loadSolver(std::stringstream* const stream) {
    int x, y;
    // reset of the stringstream
    stream->seekg(0);
    (*stream) >> x >> y;
    assert(physicalNode_.getXPos() == x && "The position does not match");
    assert(physicalNode_.getYPos() == y && "The position does not match");
    for (auto d : dirIter_) {
        (*stream) >> distributions_[d];
    }
}

void FluidSolver::resetForce() {
    force_.x = 0.0;
    force_.y = 0.0;
}

double& FluidSolver::accessDistribution(const Direction& dir) {
    return distributions_[dir];
}

void FluidSolver::rescaleDistributions(const double factor) {
    for (auto &it: this->distributions_) {
        it *= factor;
    }
}

void FluidSolver::setVelocity(Field<double> velocity) {
    velocity_ = velocity;
    assert(std::isfinite(this->velocity_.x));
    assert(std::isfinite(this->velocity_.y));
    velocityInit_ = true;
}

const Field<double>& FluidSolver::getVelocity() const {
    assert(std::isfinite(this->velocity_.x));
    assert(std::isfinite(this->velocity_.y));
    return velocity_;
}

double FluidSolver::getRho() const {
    return std::accumulate(distributions_.begin(),
            distributions_.end(), 0.0);
}

void FluidSolver::addMass(double mass) {
    this->rhoToAdd_ = 0.0;
    this->rhoToAdd_ += mass;
    LOG(UtilLib::logDEBUG4) << "mass added to the fluid. " << mass;
}

void FluidSolver::localSwap() {
    std::swap(distributions_[E], distributions_[W]);
    std::swap(distributions_[NE], distributions_[SW]);
    std::swap(distributions_[N], distributions_[S]);
    std::swap(distributions_[NW], distributions_[SE]);
}

FluidSolver::FluidSolver(const nodes::PhysicalNode& physicalNode)
    : AbstractSolver(),
      physicalNode_(physicalNode),
      distributions_(std::array<double,
                      9> {{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }
              }),
      velocityInit_(false),
      rhoToAdd_(0.0)
{}

DirectionIterator const FluidSolver::dirIter_ = DirectionIterator();
}
}  // end namespace