ChuteBottom.cc

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


#include "ChuteBottom.h"
#include "Species/LinearViscoelasticSpecies.h"
#include "Species/FrictionForceSpecies/SlidingFrictionSpecies.h"
#include "Boundaries/PeriodicBoundary.h"
#include "Walls/InfiniteWall.h"
#include <Logger.h>

ChuteBottom::ChuteBottom()
{
    constructor();
}

ChuteBottom::ChuteBottom(const DPMBase& other)
        : DPMBase(other), Chute(other)
{
    constructor();
}

ChuteBottom::ChuteBottom(const MercuryBase& other)
        : DPMBase(other), Chute(other)
{
    constructor();
}

ChuteBottom::ChuteBottom(const Mercury3D& other)
        : DPMBase(other), Chute(other)
{
    constructor();
}

ChuteBottom::ChuteBottom(const Chute& other)
        : DPMBase(other), Chute(other)
{
    constructor();
}

ChuteBottom::ChuteBottom(const ChuteBottom& other)
        : DPMBase(other), Chute(other),
          thickness_(other.thickness_), isBottomPeriodic_(other.isBottomPeriodic_)
{
    
}

void ChuteBottom::constructor()
{
    setName("roughbottom");
    fStatFile.setFileType(FileType::NO_FILE); //set to 0 for no data creation
    dataFile.setFileType(FileType::ONE_FILE);
    restartFile.setFileType(FileType::ONE_FILE);
    setThickness(2.4);
    setIsBottomPeriodic(true);
}

void ChuteBottom::makeRoughBottom(Chute& chute)
{
    // set up mini-simulation with particle inflow on the left end and horizontal chute
    // set all parameters that should be different from the original chute
    setChuteAngle(0.0);
    setInflowHeight(25. * getInflowParticleRadius());
    //~ setInflowHeight(45.*getInflowParticleRadius());
    // note: Changing the Inflow height was an attempt to make the bottom density homogeneous, but it did not have the desired effect
    setRoughBottomType(MONOLAYER_DISORDERED);
    setFixedParticleRadius(getInflowParticleRadius());

    auto species = dynamic_cast<LinearViscoelasticNormalSpecies*>(speciesHandler.getObject(0));
    if (species != nullptr)
    {
        const Mdouble mass = speciesHandler.getObject(0)->getMassFromRadius(
                0.5 * (getMinInflowParticleRadius() + getMaxInflowParticleRadius()));
        //this number is chosen to be consistent with the old implementation (with contained a bug and was replaced); originally, a factor of 10.0 (and a restitution of 0.2) was chosen.
        const Mdouble collisionTime = 13.2863 * species->getCollisionTime(mass);
        species->setCollisionTimeAndRestitutionCoefficient(collisionTime, 0.343008, mass);
        setTimeStep(0.02 * 9.41823 * collisionTime);
        logger(INFO, "Time step: %", getTimeStep());
        logger(INFO, "Species: %", *species);
    }
    else
    {
        logger(WARN, "[ChuteBottom::makeRoughBottom()] species type does not allow setting the parameters.");
    }
    
    auto species2 = dynamic_cast<SlidingFrictionSpecies*>(speciesHandler.getObject(0));
    if (species2 != nullptr)
        species2->setSlidingFrictionCoefficient(0);
    
    // set the simulation to run for 2000 time steps
    setTimeMax(getTimeStep() * 2000);
    //set_number_of_saves(2);
    setSaveCount(100);
    
    // run the simulation
    solve();

    //Find the Z-position of the highest particle in the system
    Mdouble height = 0;
    for (BaseParticle* b : particleHandler)
    {
        height = std::max(height, b->getPosition().Z);
    }
    
    // Next, all particles are removed from the system, except those with a Z-position s.t.:
    // hmax = height - maxInflowParticleRadius_;
    // hmax - thickness_ * maxInflowParticleRadius_ <  Z-position  < hmax
    // note that hmax = height - rMax
    logger(INFO, "[ChuteBottom::makeRoughBottom()] Thickness: %", thickness_);
    for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it != particleHandler.end();)
    {
        if ((*it)->getPosition().Z < height - (1.0 + thickness_) * getMaxInflowParticleRadius() ||
            (*it)->getPosition().Z > height - getMaxInflowParticleRadius())
        {
            //delete particles outside the given range
            particleHandler.removeObject((*it)->getIndex());
        }
        else
        {
            // Move remaining particles to the floor
            (*it)->move(Vec3D(0.0, 0.0, getMaxInflowParticleRadius() - height));
            
            // fix them to the bottom
            (*it)->fixParticle();
            it++;
        }
    }

    //copy the rough bottom over
    chute.particleHandler.setStorageCapacity(particleHandler.getNumberOfObjects());
    logger(INFO, "[ChuteBottom::makeRoughBottom()] Chute bottom finished, consisting of % particles",
           particleHandler.getNumberOfObjects());
    chute.particleHandler = particleHandler;
    
    for (BaseParticle* p : chute.particleHandler)
    {
        p->setSpecies(chute.speciesHandler.getObject(0));
    }

}

void ChuteBottom::setupInitialConditions()
{

    particleHandler.setStorageCapacity(static_cast<unsigned int>(std::min(3.0 * getXMax() * getYMax() * getZMax() / mathsFunc::cubic(2.0 * getInflowParticleRadius()), 1e6)));

    createBottom();

    wallHandler.clear();
    boundaryHandler.clear();
    if (isBottomPeriodic_)
    {
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin() - getInflowParticleRadius()));
        wallHandler.copyAndAddObject(w0);
        PeriodicBoundary b0;
        b0.set(Vec3D(1.0, 0.0, 0.0), getXMin(), getXMax());
        boundaryHandler.copyAndAddObject(b0);
        b0.set(Vec3D(0.0, 1.0, 0.0), getYMin(), getYMax());
        boundaryHandler.copyAndAddObject(b0);
    }
    else
    {
        InfiniteWall w0;
        w0.setSpecies(speciesHandler.getObject(0));
        w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin() - getInflowParticleRadius()));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0, 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0, getYMin(), 0));
        wallHandler.copyAndAddObject(w0);
        w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0, getYMax(), 0));
        wallHandler.copyAndAddObject(w0);
    }

    // add particles
    hGridActionsBeforeTimeLoop();
    unsigned int failed = 0;
    const unsigned int max_failed = 500;

    BaseParticle inflowParticle_;
    inflowParticle_.setSpecies(speciesHandler.getObject(0));
    inflowParticle_.setHandler(&particleHandler);
    inflowParticle_.setOrientation(Vec3D(0.0, 0.0, 0.0));
    inflowParticle_.setAngularVelocity(Vec3D(0.0, 0.0, 0.0));

    // try max_failed times to find new insertable particle
    // every time an insertion succeeds the counter failed is set to 0 again.
    while (failed <= max_failed)
    {
        inflowParticle_.setRadius(getFixedParticleRadius());
        //inflowParticle_.computeMass();

        // The position components are first stored in a Vec3D, because if you pass 
        // them directly into setPosition the compiler is allowed to change the order 
        // in which the numbers are generated
        Vec3D position;
        position.X = random.getRandomNumber(inflowParticle_.getRadius(), getXMax() - inflowParticle_.getRadius());
        position.Y = random.getRandomNumber(inflowParticle_.getRadius(), getYMax() - inflowParticle_.getRadius());
        position.Z = random.getRandomNumber(2 * inflowParticle_.getRadius(), getZMax() - inflowParticle_.getRadius());
        inflowParticle_.setPosition(position);
        inflowParticle_.setVelocity(Vec3D(0.0, 0.0, 0.0));
        
        //if the volume we want to insert the particle is still free, insert the particle
        if (checkParticleForInteraction(inflowParticle_))
        {
            particleHandler.copyAndAddObject(inflowParticle_);
            failed = 0;
        }
        else
        {
            failed++;
        }
    }
    //set_Nmax(particleHandler.getNumberOfObjects());
    logger(INFO, "[ChuteBottom::setupInitialConditions()] Number of particles created: %",
           particleHandler.getNumberOfObjects());
    
    //fix hgrid (there is still an issue when particles are polydispersed)
    //assume 1-2 levels are optimal (which is the case for mono and bidispersed) and set the cell size to min and max
    // !this is not optimal for polydispersed
    Mdouble minCell = 2. * std::min(getFixedParticleRadius(), getMinInflowParticleRadius());
    Mdouble maxCell = 2. * std::max(getFixedParticleRadius(), getMaxInflowParticleRadius());
    if ((minCell == maxCell) | (minCell == 0.))
        setHGridMaxLevels(1);
    else
        setHGridMaxLevels(2);
    //set_HGRID_cell_to_cell_ratio (1.0000000001*maxCell/minCell);
    //optimize number of buckets
    //setHGridNumberOfBucketsToPower(particleHandler.getNumberOfObjects() * 1.5);
    //end: fix hgrid

    //~ write(std::cout,false);
}

void ChuteBottom::actionsBeforeTimeStep()
{
}

Mdouble ChuteBottom::getThickness() const
{
    return thickness_;
}

void ChuteBottom::setThickness(Mdouble new_)
{
    if (new_ > 0.0)
        thickness_ = new_;
    else
    {
        logger(ERROR, "[ChuteBottom::setThickness()] thickness % negative.", new_);
    }
}

bool ChuteBottom::getIsBottomPeriodic() const
{
    return isBottomPeriodic_;
}

void ChuteBottom::setIsBottomPeriodic(bool isBottomPeriodic)
{
    isBottomPeriodic_ = isBottomPeriodic;
}