diff --git a/src/lbm/SuViscoelasticity.h b/src/lbm/SuViscoelasticity.h
index 481d9151e8dca00e9217ef75441511bb4954bd97..15e1a7195fd036bc559c1bd9f381574aa5cf666b 100644
--- a/src/lbm/SuViscoelasticity.h
+++ b/src/lbm/SuViscoelasticity.h
@@ -16,27 +16,17 @@
//! \file SuViscoelasticity.h
//! \ingroup lbm
//! \author Cameron Stewart <cstewart@icp.uni-stuttgart.de>
-//! \brief Oldroyd-B viscoelasticity modeled my Su et al. Phys. Rev. E, 2013
+//! \brief Oldroyd-B viscoelasticity extended from Su et al. Phys. Rev. E, 2013
//
//======================================================================================================================
-#include "field/GhostLayerField.h"
-#include "../field/AddToStorage.h"
-#include "domain_decomposition/IBlock.h"
+
+
+#include "blockforest/communication/UniformBufferedScheme.h"
#include "core/math/Matrix3.h"
-#include "core/DataTypes.h"
-#include "stencil/Directions.h"
-#include "stencil/D2Q4.h"
+#include <field/AddToStorage.h>
+#include "field/GhostLayerField.h"
#include "field/communication/PackInfo.h"
#include "lbm/field/PdfField.h"
-#include "blockforest/communication/UniformBufferedScheme.h"
-#include "lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h"
-
-#include "stencil/D2Q9.h"
-#include "stencil/D3Q6.h"
-
-#include <iostream>
-#include <math.h>
-#include <field/all.h>
namespace walberla {
namespace lbm {
@@ -101,8 +91,9 @@ public:
Matrix3<real_t> relstr = Matrix3<real_t>(0.0);
bool nearBoundaryFlag = false;
- // If cell is a fluid cell then calculate the stress tensor
+ // if cell is a fluid cell then calculate the stress tensor
if (boundaryHandling->isDomain(cell)) {
+ // check if near a wall
if (boundaryHandling->isNearBoundary(cell)) {
nearBoundaryFlag = true;
}
@@ -119,7 +110,7 @@ public:
Cell cell2 = cell1 - *d;
Cell cell3 = cell + *d;
- // Check if using compressible of incompressible algorithm
+ // check if using compressible of incompressible algorithm
if (compressibleFlag_) {
if (nearBoundaryFlag) {
if (boundaryHandling->isDomain(cell1)) {
@@ -151,9 +142,9 @@ public:
}
}
}
- // Compute velocity gradient
Matrix3<real_t> gradu = Matrix3<real_t>(0.0);
+ // compute velocity gradient
for (auto d = LatticeModel_T::Stencil::beginNoCenter(); d.direction() != stencil::NW; ++d) {
for (uint_t a = 0; a < LatticeModel_T::Stencil::D; ++a) {
for (uint_t b = 0; b < LatticeModel_T::Stencil::D; ++b) {
@@ -169,14 +160,14 @@ public:
}
auto graduT = gradu.getTranspose();
- // Equation 16 from Su 2013
+ // equation 16 from Su 2013
relstr = stressOld->get(cell) * gradu + graduT * stressOld->get(cell);
if(eta_p_ > real_t(0)) {
relstr += ((gradu + graduT) * eta_p_ - stressOld->get(cell)) * inv_lambda_p_;
}
- // Equation 23 from Su 2013
+ // equation 23 from Su 2013
stressNew->get(cell) = stressOld->get(cell) + (stress1 * real_c(2.0) - stress2 * real_c(1.5) - stress3 * real_c(0.5)) * delta_t_ * (real_c(1) / pdf->getDensity(cell)) +
relstr * delta_t_;
}
@@ -198,8 +189,8 @@ public:
WALBERLA_ASSERT_GREATER_EQUAL( velocity->nrOfGhostLayers(), 1 );
- // Update velocity field for all fluid cells
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(velocity,{
+ // update velocity field for all fluid cells
+ WALBERLA_FOR_ALL_CELLS_XYZ(velocity,{
Cell cell(x, y, z);
if( boundaryHandling->isDomain(cell) ) {
velocity->get( cell ) = pdf->getVelocity(x, y, z);
diff --git a/tests/lbm/Su.prm b/tests/lbm/Su.prm
index 6c215a2cb4f92646f3e2030abe09d398246fbea4..88a9bfe3baaa384708a746f6ac8a98dd8aa42727 100644
--- a/tests/lbm/Su.prm
+++ b/tests/lbm/Su.prm
@@ -8,6 +8,7 @@ Parameters
period 1;
L 11;
H 33;
+ blockSize 11;
uMax 1.28737;
t0 370;
t1 326;
diff --git a/tests/lbm/SuViscoelasticityTest.cpp b/tests/lbm/SuViscoelasticityTest.cpp
index cecac2e28b68c12455de2e398a8ed8d93b0c85f9..fa786bcc25a2c81b804f6849290b8f63b5e17172 100644
--- a/tests/lbm/SuViscoelasticityTest.cpp
+++ b/tests/lbm/SuViscoelasticityTest.cpp
@@ -16,49 +16,34 @@
//! \file SuViscoelasticityTest.cpp
//! \ingroup lbm
//! \author Cameron Stewart <cstewart@icp.uni-stuttgart.de>
-//! \brief D2Q9 Poisuielle flow simulation with viscoelasticity
+//! \brief D2Q9 Poiseuille flow simulation with Oldroyd-B viscoelasticity - Checks against analytical formula and
+//! reference data
//
//
//======================================================================================================================
-
#include "blockforest/Initialization.h"
#include "blockforest/communication/UniformBufferedScheme.h"
#include "boundary/BoundaryHandling.h"
-#include "core/Abort.h"
#include "core/Environment.h"
-#include "core/debug/TestSubsystem.h"
-#include "core/logging/Logging.h"
-#include "core/timing/RemainingTimeLogger.h"
#include "core/SharedFunctor.h"
-#include "core/Filesystem.h"
#include "domain_decomposition/SharedSweep.h"
#include "field/communication/PackInfo.h"
-#include "field/StabilityChecker.h"
-
-#include "geometry/initializer/BoundaryFromDomainBorder.h"
-#include "geometry/InitBoundaryHandling.h"
-#include "lbm/lattice_model/D2Q9.h"
-#include "lbm/field/AddToStorage.h"
#include "lbm/boundary/all.h"
-#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/lattice_model/D2Q9.h"
#include "lbm/lattice_model/ForceModel.h"
#include "lbm/SuViscoelasticity.h"
#include "lbm/sweeps/CellwiseSweep.h"
-#include "lbm/vtk/all.h"
#include "timeloop/SweepTimeloop.h"
-#include <math.h>
#include <iostream>
-#include <iomanip>
-#include <vector>
-#include <lbm/boundary/factories/ExtendedBoundaryHandlingFactory.h>
using namespace walberla;
@@ -66,7 +51,6 @@ using namespace walberla;
// TYPEDEFS //
//////////////
-
typedef GhostLayerField< Vector3<real_t>, 1> ForceField_T;
typedef lbm::collision_model::TRT CollisionModel_T;
@@ -261,21 +245,21 @@ int main(int argc, char ** argv ){
auto parameters = env.config()->getOneBlock( "Parameters" );
// extract some constants from the parameters
- const real_t eta_s = parameters.getParameter< real_t > ( "eta_s", real_c( 1.4 ) );
- const real_t force = parameters.getParameter< real_t > ( "force", real_c( 1 ) );
- const real_t eta_p = parameters.getParameter< real_t > ( "eta_p", real_c( 1.0 ) );
- const real_t lambda_p = parameters.getParameter< real_t > ( "lambda_p", real_c( 10 ) );
- const uint_t period = parameters.getParameter< uint_t > ( "period", uint_c( 1 ) );
- const real_t L = parameters.getParameter< real_t > ( "L", real_c( 10 ) );
- const real_t H = parameters.getParameter< real_t > ( "H", real_c( 30 ) );
- const uint_t blockSize = parameters.getParameter< uint_t > ( "blockSize", uint_c( 11 ) );
- const uint_t timesteps = parameters.getParameter< uint_t > ( "timesteps", uint_c( 10 ) );
+ const real_t eta_s = parameters.getParameter< real_t > ("eta_s");
+ const real_t force = parameters.getParameter< real_t > ("force");
+ const real_t eta_p = parameters.getParameter< real_t > ("eta_p");
+ const real_t lambda_p = parameters.getParameter< real_t > ("lambda_p");
+ const uint_t period = parameters.getParameter< uint_t > ("period");
+ const real_t L = parameters.getParameter< real_t > ("L");
+ const real_t H = parameters.getParameter< real_t > ("H");
+ const uint_t blockSize = parameters.getParameter< uint_t > ("blockSize");
+ const uint_t timesteps = parameters.getParameter< uint_t > ("timesteps");
// reference data
- const real_t uExpected = force*H*H/(8.0*(eta_s + eta_p));
- const real_t uMax = parameters.getParameter< real_t > ("uMax");
- const real_t t0 = parameters.getParameter< real_t > ("t0");
- const real_t t1 = parameters.getParameter< real_t > ("t1");
+ const real_t uExpected = force*H*H/(8.0*(eta_s + eta_p));
+ const real_t uMax = parameters.getParameter< real_t > ("uMax");
+ const real_t t0 = parameters.getParameter< real_t > ("t0");
+ const real_t t1 = parameters.getParameter< real_t > ("t1");
// create fields
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field", FieldGhostLayers);
@@ -297,6 +281,7 @@ int main(int argc, char ** argv ){
communication.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldId ) );
auto testData = make_shared< TestData >(TestData(timeloop, blocks, pdfFieldId, stressId, timesteps, blockSize, L, H, uExpected));
+ // structure timeloop
timeloop.add() << BeforeFunction( communication, "communication" )
<< Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingId ), "boundary handling" );
timeloop.add() << BeforeFunction( lbm::viscoelastic::Su<LatticeModel_T, BoundaryHandling_T>(blocks, forceFieldId, pdfFieldId, boundaryHandlingId, stressId, stressOldId, velocityId,
@@ -315,18 +300,18 @@ int main(int argc, char ** argv ){
real_t errt0 = abs(testData->getT0() - t0)/t0;
real_t errt1 = abs(testData->getT1() - t1)/t1;
- std::cout << "Relative Errors" << std::endl << std::endl;
- std::cout << "Steady State Velocity: " << errSteady << std::endl;
- std::cout << "Maximum Velocity: " << errMax << std::endl;
- std::cout << "Time of Maximum: " << errt0 << std::endl;
- std::cout << "Decay Time: " << errt1 << std::endl << std::endl;
+ WALBERLA_LOG_RESULT("Steady State Velocity Error: " << errSteady );
+ WALBERLA_LOG_RESULT("Maximum Velocity Error: " << errMax );
+ WALBERLA_LOG_RESULT("Time of Maximum Error: " << errt0 );
+ WALBERLA_LOG_RESULT("Decay Time Error: " << errt1 );
+ // check that errors < 1%
if (errSteady < 0.01 && errMax < 0.01 && errt0 < 0.01 && errt1 < 0.01){
- std::cout << "Success" << std::endl;
+ WALBERLA_LOG_RESULT("Success" );
return EXIT_SUCCESS;
}
else {
- std::cout << "Failure" << std::endl;
+ WALBERLA_LOG_RESULT("Failure" );
return EXIT_FAILURE;
}