Skip to content
Snippets Groups Projects
Commit 3fa0d2e2 authored by Lukas Werner's avatar Lukas Werner
Browse files

#9 deleted no longer necessary files

parent 0fb9318a
No related merge requests found
waLBerla_add_executable ( NAME MESA_PD_VIRTUALMASS
FILES VirtualMass.cpp
DEPENDS mesa_pd core )
waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_VIRTUALMASS )
waLBerla_add_executable ( NAME MESA_PD_APP_UNSETTLINGSPHERE_VM
FILES UnsettlingSphere.cpp
DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field vtk )
waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_UNSETTLINGSPHERE_VM )
waLBerla_add_executable ( NAME MESA_PD_APP_FIXEDROTATINGSPHERE_VM
FILES FixedRotatingSphere.cpp
DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field vtk )
waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_FIXEDROTATINGSPHERE_VM )
waLBerla_add_executable ( NAME MESA_PD_APP_SMALLUNSETTLINGSPHERE_VM
FILES SmallUnsettlingSphere.cpp
DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field vtk )
waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_SMALLUNSETTLINGSPHERE_VM )
DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field vtk )
\ No newline at end of file
This diff is collapsed.
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file VirtualMass.cpp
//! \author Lukas Werner <lks.werner@fau.de>
//
//======================================================================================================================
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
#include <mesa_pd/kernel/AddVirtualForceAndTorque.h>
#include <mesa_pd/kernel/VirtualMass.h>
#include <core/Environment.h>
#include <core/logging/Logging.h>
#include <iostream>
namespace walberla {
using namespace walberla::mesa_pd;
class ParticleAccessorWithShape : public data::ParticleAccessor
{
public:
ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
: ParticleAccessor(ps)
, ss_(ss)
{}
const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass() = v;}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
void setInvInertiaBF(const size_t p_idx, const Mat3& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF() = v;}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
walberla::real_t getVolume(const size_t p_idx) const { return getShape(p_idx)->getVolume(); }
private:
std::shared_ptr<data::ShapeStorage> ss_;
};
/**
* Accessor which takes the virtual mass for calls to invMass and getInvInertiaBF into account.
*/
class ParticleAccessorWithShapeAndVirtualMass : public ParticleAccessorWithShape {
public:
ParticleAccessorWithShapeAndVirtualMass(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
: ParticleAccessorWithShape(ps, ss)
{}
/**
* Returns the inverse mass of the body by also calculating and including its virtual mass.
*
* \param p_idx index of particle
*/
walberla::real_t getInvMass(const size_t p_idx) const {
return getInvMassIncludingVirtual(p_idx);
}
const auto getInvInertiaBF(const size_t p_idx) const {
return getInvInertiaBFIncludingVirtual(p_idx);
}
walberla::real_t& getInvMassRef(const size_t p_idx) {
WALBERLA_UNUSED(p_idx);
WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
return ParticleAccessorWithShape::getInvMassRef(p_idx);
}
void setInvMass(const size_t p_idx, const walberla::real_t& v) {
WALBERLA_UNUSED(p_idx);
WALBERLA_UNUSED(v);
WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
}
auto& getInvInertiaBFRef(const size_t p_idx) {
WALBERLA_UNUSED(p_idx);
WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
return ParticleAccessorWithShape::getInvInertiaBFRef(p_idx);
}
void setInvInertiaBF(const size_t p_idx, const Mat3& v) {
WALBERLA_UNUSED(p_idx);
WALBERLA_UNUSED(v);
WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
}
};
int main( int argc, char ** argv )
{
Environment env(argc, argv);
WALBERLA_UNUSED(env);
mpi::MPIManager::instance()->useWorldComm();
auto ps = std::make_shared<data::ParticleStorage>(2);
auto ss = std::make_shared<data::ShapeStorage>();
ParticleAccessorWithShape accessor(ps, ss);
// accessor for the virtual mass approach
ParticleAccessorWithShapeAndVirtualMass virtualMassAccessor(ps, ss);
//initialize particles
const auto radius = real_t(0.6);
//const auto spacing = real_t(1.0);
auto smallSphere = ss->create<data::Sphere>( radius );
ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
ps->create();
ps->create();
auto force = Vec3(1,2,3);
auto torque = Vec3(1,2,3);
ps->forEachParticle(false, [force, torque](size_t idx, auto& ac) {
const auto linVel = Vec3(1,2,3);
const auto angVel = Vec3(1,2,3);
const auto invMass = real_t(1.23456);
const auto invInertiaBF = Mat3(
real_t(1.23456), real_t(0), real_t(0),
real_t(0), real_t(1.23456), real_t(0),
real_t(0), real_t(0), real_t(1.23456));
ac.setLinearVelocity(idx, linVel);
ac.setAngularVelocity(idx, angVel);
ac.setForce(idx, force);
ac.setTorque(idx, torque);
ac.setInvMass(idx, invMass);
ac.setInvInertiaBF(idx, invInertiaBF);
}, accessor);
// indices of test particles
const size_t plainIdx = 0; // particle using the default integration
const size_t virtualIdx = 1; // particle using the virtual mass approach
//init kernels
const auto dt = real_t(1);
// constants required for virtual mass approach
const auto C_v = real_t(0.5);
const auto C_v_omega = C_v;
const auto fluidDensity = real_t(1);
//kernel::ExplicitEuler integrator( dt );
kernel::ExplicitEulerWithShape integrator( dt );
kernel::VirtualMass virtualMass(C_v, C_v_omega, fluidDensity);
kernel::AddVirtualForceAndTorque addVirtualForceAndTorque;
ps->forEachParticle(false, [virtualMass](size_t idx, auto& ac) {
virtualMass(idx, ac);
}, accessor);
for (int i = 0; i < 3; ++i) {
// integration of particle with no virtual mass added
accessor.setForce(plainIdx, force);
accessor.setTorque(plainIdx, torque);
integrator(plainIdx, accessor);
// integration of particle with virtual mass added
accessor.setForce(virtualIdx, force);
accessor.setTorque(virtualIdx, torque);
addVirtualForceAndTorque(virtualIdx, accessor);
integrator(virtualIdx, virtualMassAccessor);
force += Vec3(0.1, 0.1, 0.1);
torque += Vec3(0.1, 0.1, 0.1);
WALBERLA_LOG_INFO("p vel: " << accessor.getLinearVelocity(plainIdx));
WALBERLA_LOG_INFO("p angvel: " << accessor.getAngularVelocity(plainIdx));
WALBERLA_LOG_INFO("v vel: " << accessor.getLinearVelocity(virtualIdx));
WALBERLA_LOG_INFO("v angvel: " << accessor.getAngularVelocity(virtualIdx));
WALBERLA_ASSERT_FLOAT_EQUAL(accessor.getLinearVelocity(plainIdx), accessor.getLinearVelocity(virtualIdx));
WALBERLA_ASSERT_FLOAT_EQUAL(accessor.getAngularVelocity(plainIdx), accessor.getAngularVelocity(virtualIdx));
}
return EXIT_SUCCESS;
}
} //namespace walberla
int main( int argc, char ** argv )
{
return walberla::main(argc, argv);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment