diff --git a/CMake/cmake_config_files/AMATERASU.config.cmake b/CMake/cmake_config_files/AMATERASU.config.cmake
index 4cf53d794e743b61aa18a0565f97f6eb5fbf723b..e65b197cba1a5d70007144c09cd80c85c33a0696 100644
--- a/CMake/cmake_config_files/AMATERASU.config.cmake
+++ b/CMake/cmake_config_files/AMATERASU.config.cmake
@@ -9,6 +9,6 @@
#################################################################################
# METIS
#################################################################################
-SET(METIS_INCLUDEDIR "C:/Users/geier/Documents/metis-5.1.0/include")
-SET(METIS_DEBUG_LIBRARY "C:/Users/geier/Documents/metis-5.1.0/build/libmetis/Debug/metis.lib")
-SET(METIS_RELEASE_LIBRARY "C:/Users/geier/Documents/metis-5.1.0/build/libmetis/Release/metis.lib")
+#SET(METIS_INCLUDEDIR "C:/Users/geier/Documents/metis-5.1.0/include")
+#SET(METIS_DEBUG_LIBRARY "C:/Users/geier/Documents/metis-5.1.0/build/libmetis/Debug/metis.lib")
+#SET(METIS_RELEASE_LIBRARY "C:/Users/geier/Documents/metis-5.1.0/build/libmetis/Release/metis.lib")
diff --git a/CMake/cmake_config_files/BOMBADIL.config.cmake b/CMake/cmake_config_files/BOMBADIL.config.cmake
index 9a205596d22cfac18baef3131871b9f73e1b4f9c..c0c5cf2f08b2cc925be248441c99a336160fd1bd 100644
--- a/CMake/cmake_config_files/BOMBADIL.config.cmake
+++ b/CMake/cmake_config_files/BOMBADIL.config.cmake
@@ -29,15 +29,15 @@ set(VTK_DIR "d:/Tools/VTK/build/VTK-8.0.0")
#################################################################################
# METIS
#################################################################################
-IF(${USE_METIS})
- SET(METIS_INCLUDEDIR "d:/Tools/metis-5.1.0/include")
- SET(METIS_DEBUG_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Debug/metis.lib")
- SET(METIS_RELEASE_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Release/metis.lib")
+#IF(${USE_METIS})
+# SET(METIS_INCLUDEDIR "d:/Tools/metis-5.1.0/include")
+# SET(METIS_DEBUG_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Debug/metis.lib")
+# SET(METIS_RELEASE_LIBRARY "d:/Tools/metis-5.1.0/build/libmetis/Release/metis.lib")
# SET(METIS_INCLUDEDIR "/mnt/d/Tools/metis-5.1.0/include")
# SET(METIS_DEBUG_LIBRARY "/mnt/d/Tools/metis-5.1.0/build/Linux-x86_64/libmetis/libmetis.a")
# SET(METIS_RELEASE_LIBRARY "/mnt/d/Tools/metis-5.1.0/build/Linux-x86_64/libmetis/libmetis.a")
-ENDIF()
+#ENDIF()
#################################################################################
# PE
diff --git a/apps/cpu/Multiphase/Multiphase.cfg b/apps/cpu/Multiphase/Multiphase.cfg
index 0d3e5633ea5b4c718eadd5fba04b89abc38497ed..34ef25bbc20475c16e011191f652ed1a1de3a2f4 100644
--- a/apps/cpu/Multiphase/Multiphase.cfg
+++ b/apps/cpu/Multiphase/Multiphase.cfg
@@ -44,5 +44,5 @@ restartStep = 100000
cpStart = 100000
cpStep = 100000
-outTime = 1
-endTime = 200000000
\ No newline at end of file
+outTime = 100
+endTime = 2000
\ No newline at end of file
diff --git a/apps/cpu/Multiphase/Multiphase.cpp b/apps/cpu/Multiphase/Multiphase.cpp
index 1973ae061c02118c3a08210ece90f1b6905ef81a..f122e90157d8ecbdc7128a156edba407414ad2d0 100644
--- a/apps/cpu/Multiphase/Multiphase.cpp
+++ b/apps/cpu/Multiphase/Multiphase.cpp
@@ -81,7 +81,7 @@ void run(string configname)
SPtr<LBMKernel> kernel;
- kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
+ kernel = SPtr<LBMKernel>(new MultiphaseScratchCumulantLBMKernel());
kernel->setWithForcing(true);
kernel->setForcingX1(0.0);
@@ -471,7 +471,7 @@ void run(string configname)
grid->accept(setConnsVisitor);
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
- SPtr<WriteMacroscopicQuantitiesCoProcessor> pp(new WriteMacroscopicQuantitiesCoProcessor(
+ SPtr<WriteMultiphaseQuantitiesCoProcessor> pp(new WriteMultiphaseQuantitiesCoProcessor(
grid, visSch, pathname, WbWriterVtkXmlASCII::getInstance(), conv, comm));
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
diff --git a/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index d60189b0747541e5e02fb38752df4581ee754b33..564cc2768469d0fdb83647126d2c94b6314d146b 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -204,6 +204,7 @@
#include <CoProcessors/MPIIORestartCoProcessor.h>
#include <CoProcessors/MicrophoneArrayCoProcessor.h>
#include <WriteThixotropyQuantitiesCoProcessor.h>
+#include <WriteMultiphaseQuantitiesCoProcessor.h>
#include <IntegrateValuesHelper.h>
//#include <LBM/D3Q27CompactInterpolationProcessor.h>
@@ -239,6 +240,7 @@
#include <LBM/RheologyK17LBMKernel.h>
#include <LBM/RheologyPowellEyringModelLBMKernel.h>
#include <LBM/MultiphaseCumulantLBMKernel.h>
+#include <LBM/MultiphaseScratchCumulantLBMKernel.h>
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6c418c11ad3cbddcf13c64f93ffdc156aa766912
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.cpp
@@ -0,0 +1,426 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids 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.
+//
+// VirtualFluids 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 VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file WriteMultiphaseQuantitiesCoProcessor.cpp
+//! \ingroup CoProcessors
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#include "WriteMultiphaseQuantitiesCoProcessor.h"
+#include "BCProcessor.h"
+#include "LBMKernel.h"
+#include <string>
+#include <vector>
+
+#include "BCArray3D.h"
+#include "Block3D.h"
+#include "Communicator.h"
+#include "DataSet3D.h"
+#include "Grid3D.h"
+#include "LBMUnitConverter.h"
+#include "UbScheduler.h"
+#include "basics/writer/WbWriterVtkXmlASCII.h"
+
+WriteMultiphaseQuantitiesCoProcessor::WriteMultiphaseQuantitiesCoProcessor() = default;
+//////////////////////////////////////////////////////////////////////////
+WriteMultiphaseQuantitiesCoProcessor::WriteMultiphaseQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
+ const std::string &path,
+ WbWriter *const writer,
+ SPtr<LBMUnitConverter> conv,
+ SPtr<Communicator> comm)
+ : CoProcessor(grid, s), path(path), writer(writer), conv(conv), comm(comm)
+{
+ gridRank = comm->getProcessID();
+ minInitLevel = this->grid->getCoarsestInitializedLevel();
+ maxInitLevel = this->grid->getFinestInitializedLevel();
+
+ blockVector.resize(maxInitLevel + 1);
+
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
+ grid->getBlocks(level, gridRank, true, blockVector[level]);
+ }
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::init()
+{}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::process(double step)
+{
+ if (scheduler->isDue(step))
+ collectData(step);
+
+ UBLOG(logDEBUG3, "WriteMultiphaseQuantitiesCoProcessor::update:" << step);
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::collectData(double step)
+{
+ int istep = static_cast<int>(step);
+
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
+ for (SPtr<Block3D> block : blockVector[level])
+ {
+ if (block)
+ {
+ addDataMQ(block);
+ }
+ }
+ }
+
+ std::string pfilePath, partPath, subfolder, cfilePath;
+
+ subfolder = "mq" + UbSystem::toString(istep);
+ pfilePath = path + "/mq/" + subfolder;
+ cfilePath = path + "/mq/mq_collection";
+ partPath = pfilePath + "/mq" + UbSystem::toString(gridRank) + "_" + UbSystem::toString(istep);
+
+ std::string partName = writer->writeOctsWithNodeData(partPath, nodes, cells, datanames, data);
+ size_t found = partName.find_last_of("/");
+ std::string piece = partName.substr(found + 1);
+ piece = subfolder + "/" + piece;
+
+ std::vector<std::string> cellDataNames;
+ std::vector<std::string> pieces = comm->gather(piece);
+ if (comm->getProcessID() == comm->getRoot()) {
+ std::string pname =
+ WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
+ found = pname.find_last_of("/");
+ piece = pname.substr(found + 1);
+
+ std::vector<std::string> filenames;
+ filenames.push_back(piece);
+ if (step == CoProcessor::scheduler->getMinBegin())
+ {
+ WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
+ } else
+ {
+ WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
+ }
+ UBLOG(logINFO, "WriteMultiphaseQuantitiesCoProcessor step: " << istep);
+ }
+
+ clearData();
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::clearData()
+{
+ nodes.clear();
+ cells.clear();
+ datanames.clear();
+ data.clear();
+}
+
+//////////////////////////////////////////////////////////////////////////
+void WriteMultiphaseQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
+{
+ using namespace D3Q27System;
+ using namespace UbMath;
+
+ double level = (double)block->getLevel();
+
+ // Diese Daten werden geschrieben:
+ datanames.resize(0);
+ datanames.push_back("Phi");
+ datanames.push_back("Vx");
+ datanames.push_back("Vy");
+ datanames.push_back("Vz");
+ datanames.push_back("P1");
+
+ data.resize(datanames.size());
+
+ SPtr<LBMKernel> kernel = dynamicPointerCast<LBMKernel>(block->getKernel());
+ SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+ SPtr<DistributionArray3D> distributionsF = kernel->getDataSet()->getFdistributions();
+ SPtr<DistributionArray3D> distributionsH = kernel->getDataSet()->getHdistributions();
+ SPtr<PhaseFieldArray3D> divU = kernel->getDataSet()->getPhaseField();
+
+ LBMReal f[D3Q27System::ENDF + 1];
+ LBMReal phi[D3Q27System::ENDF + 1];
+ LBMReal vx1, vx2, vx3, rho, p1, beta, kappa;
+ LBMReal densityRatio = kernel->getDensityRatio();
+
+ kernel->getMultiphaseModelParameters(beta, kappa);
+ LBMReal phiL = kernel->getPhiL();
+ LBMReal phiH = kernel->getPhiH();
+
+ // knotennummerierung faengt immer bei 0 an!
+ int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+
+ if (block->getKernel()->getCompressible()) {
+ calcMacros = &D3Q27System::calcCompMacroscopicValues;
+ } else {
+ calcMacros = &D3Q27System::calcIncompMacroscopicValues;
+ }
+
+ // int minX1 = 0;
+ // int minX2 = 0;
+ // int minX3 = 0;
+
+ int maxX1 = (int)(distributionsF->getNX1());
+ int maxX2 = (int)(distributionsF->getNX2());
+ int maxX3 = (int)(distributionsF->getNX3());
+
+ int minX1 = 0;
+ int minX2 = 0;
+ int minX3 = 0;
+
+ // int maxX1 = (int)(distributions->getNX1());
+ // int maxX2 = (int)(distributions->getNX2());
+ // int maxX3 = (int)(distributions->getNX3());
+
+ // nummern vergeben und node vector erstellen + daten sammeln
+ CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr phaseField1(
+ new CbArray3D<LBMReal, IndexerX3X2X1>(maxX1, maxX2, maxX3, -999.0));
+
+ for (size_t ix3 = minX3; ix3 < maxX3; ix3++) {
+ for (size_t ix2 = minX2; ix2 < maxX2; ix2++) {
+ for (size_t ix1 = minX1; ix1 < maxX1; ix1++) {
+ if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3)) {
+ distributionsH->getDistribution(f, ix1, ix2, ix3);
+ (*phaseField1)(ix1, ix2, ix3) =
+ ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
+ (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
+ ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
+ ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+ }
+ }
+ }
+ }
+
+ maxX1 -= 2;
+ maxX2 -= 2;
+ maxX3 -= 2;
+
+ // maxX1 -= 1;
+ // maxX2 -= 1;
+ // maxX3 -= 1;
+
+ // D3Q27BoundaryConditionPtr bcPtr;
+ int nr = (int)nodes.size();
+ LBMReal dX1_phi;
+ LBMReal dX2_phi;
+ LBMReal dX3_phi;
+ LBMReal mu;
+
+ for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
+ for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
+ for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
+ if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3)) {
+ int index = 0;
+ nodeNumbers(ix1, ix2, ix3) = nr++;
+ Vector3D worldCoordinates = grid->getNodeCoordinates(block, ix1, ix2, ix3);
+ nodes.push_back(UbTupleFloat3(float(worldCoordinates[0]), float(worldCoordinates[1]),
+ float(worldCoordinates[2])));
+
+ phi[REST] = (*phaseField1)(ix1, ix2, ix3);
+
+ if ((ix1 == 0) || (ix2 == 0) || (ix3 == 0)) {
+ dX1_phi = 0.0;
+ dX2_phi = 0.0;
+ dX3_phi = 0.0;
+ mu = 0.0;
+ // vx1 = 0.0;
+ // vx2 = 0.0;
+ // vx3 = 0.0;
+ } else {
+ phi[E] = (*phaseField1)(ix1 + DX1[E], ix2 + DX2[E], ix3 + DX3[E]);
+ phi[N] = (*phaseField1)(ix1 + DX1[N], ix2 + DX2[N], ix3 + DX3[N]);
+ phi[T] = (*phaseField1)(ix1 + DX1[T], ix2 + DX2[T], ix3 + DX3[T]);
+ phi[W] = (*phaseField1)(ix1 + DX1[W], ix2 + DX2[W], ix3 + DX3[W]);
+ phi[S] = (*phaseField1)(ix1 + DX1[S], ix2 + DX2[S], ix3 + DX3[S]);
+ phi[B] = (*phaseField1)(ix1 + DX1[B], ix2 + DX2[B], ix3 + DX3[B]);
+ phi[NE] = (*phaseField1)(ix1 + DX1[NE], ix2 + DX2[NE], ix3 + DX3[NE]);
+ phi[NW] = (*phaseField1)(ix1 + DX1[NW], ix2 + DX2[NW], ix3 + DX3[NW]);
+ phi[TE] = (*phaseField1)(ix1 + DX1[TE], ix2 + DX2[TE], ix3 + DX3[TE]);
+ phi[TW] = (*phaseField1)(ix1 + DX1[TW], ix2 + DX2[TW], ix3 + DX3[TW]);
+ phi[TN] = (*phaseField1)(ix1 + DX1[TN], ix2 + DX2[TN], ix3 + DX3[TN]);
+ phi[TS] = (*phaseField1)(ix1 + DX1[TS], ix2 + DX2[TS], ix3 + DX3[TS]);
+ phi[SW] = (*phaseField1)(ix1 + DX1[SW], ix2 + DX2[SW], ix3 + DX3[SW]);
+ phi[SE] = (*phaseField1)(ix1 + DX1[SE], ix2 + DX2[SE], ix3 + DX3[SE]);
+ phi[BW] = (*phaseField1)(ix1 + DX1[BW], ix2 + DX2[BW], ix3 + DX3[BW]);
+ phi[BE] = (*phaseField1)(ix1 + DX1[BE], ix2 + DX2[BE], ix3 + DX3[BE]);
+ phi[BS] = (*phaseField1)(ix1 + DX1[BS], ix2 + DX2[BS], ix3 + DX3[BS]);
+ phi[BN] = (*phaseField1)(ix1 + DX1[BN], ix2 + DX2[BN], ix3 + DX3[BN]);
+ phi[BSW] = (*phaseField1)(ix1 + DX1[BSW], ix2 + DX2[BSW], ix3 + DX3[BSW]);
+ phi[BSE] = (*phaseField1)(ix1 + DX1[BSE], ix2 + DX2[BSE], ix3 + DX3[BSE]);
+ phi[BNW] = (*phaseField1)(ix1 + DX1[BNW], ix2 + DX2[BNW], ix3 + DX3[BNW]);
+ phi[BNE] = (*phaseField1)(ix1 + DX1[BNE], ix2 + DX2[BNE], ix3 + DX3[BNE]);
+ phi[TNE] = (*phaseField1)(ix1 + DX1[TNE], ix2 + DX2[TNE], ix3 + DX3[TNE]);
+ phi[TNW] = (*phaseField1)(ix1 + DX1[TNW], ix2 + DX2[TNW], ix3 + DX3[TNW]);
+ phi[TSE] = (*phaseField1)(ix1 + DX1[TSE], ix2 + DX2[TSE], ix3 + DX3[TSE]);
+ phi[TSW] = (*phaseField1)(ix1 + DX1[TSW], ix2 + DX2[TSW], ix3 + DX3[TSW]);
+ dX1_phi = 0.0 * gradX1_phi(phi);
+ dX2_phi = 0.0 * gradX2_phi(phi);
+ dX3_phi = 0.0 * gradX3_phi(phi);
+ mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi(phi);
+ }
+
+ distributionsF->getDistribution(f, ix1, ix2, ix3);
+ LBMReal dU = (*divU)(ix1, ix2, ix3);
+
+ LBMReal rhoH = 1.0;
+ LBMReal rhoL = 1.0 / densityRatio;
+ // LBMReal rhoToPhi = (1.0 - 1.0/densityRatio);
+ LBMReal rhoToPhi = (rhoH - rhoL) / (phiH - phiL);
+
+ // rho = phi[ZERO] + (1.0 - phi[ZERO])*1.0/densityRatio;
+ rho = rhoH + rhoToPhi * (phi[REST] - phiH);
+
+ vx1 =
+ ((((f[TNE] - f[BSW]) + (f[TSE] - f[BNW])) + ((f[BSE] - f[TNW]) + (f[BNE] - f[TSW]))) +
+ (((f[BE] - f[TW]) + (f[TE] - f[BW])) + ((f[SE] - f[NW]) + (f[NE] - f[SW]))) + (f[E] - f[W])) /
+ (rho * c1o3) +
+ mu * dX1_phi / (2 * rho);
+
+ vx2 =
+ ((((f[TNE] - f[BSW]) + (f[BNW] - f[TSE])) + ((f[TNW] - f[BSE]) + (f[BNE] - f[TSW]))) +
+ (((f[BN] - f[TS]) + (f[TN] - f[BS])) + ((f[NW] - f[SE]) + (f[NE] - f[SW]))) + (f[N] - f[S])) /
+ (rho * c1o3) +
+ mu * dX2_phi / (2 * rho);
+
+ vx3 =
+ ((((f[TNE] - f[BSW]) + (f[TSE] - f[BNW])) + ((f[TNW] - f[BSE]) + (f[TSW] - f[BNE]))) +
+ (((f[TS] - f[BN]) + (f[TN] - f[BS])) + ((f[TW] - f[BE]) + (f[TE] - f[BW]))) + (f[T] - f[B])) /
+ (rho * c1o3) +
+ mu * dX3_phi / (2 * rho);
+
+ p1 = (((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
+ (((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
+ ((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
+ ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST]) +
+ (vx1 * rhoToPhi * dX1_phi * c1o3 + vx2 * rhoToPhi * dX2_phi * c1o3 +
+ vx3 * rhoToPhi * dX3_phi * c1o3) /
+ 2.0;
+ p1 = rho * p1 * c1o3;
+
+ // calcMacros(f,rho,vx1,vx2,vx3);
+ // tempfield(ix1,ix2,ix3)= ; Added by HS
+ // double press = D3Q27System::calcPress(f,rho,vx1,vx2,vx3);
+
+ if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1))
+ UB_THROW(UbException(
+ UB_EXARGS, "vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+ block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+ UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+ // vx1=999.0;
+ if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2))
+ UB_THROW(UbException(
+ UB_EXARGS, "vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+ block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+ UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+ // vx2=999.0;
+ if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3))
+ UB_THROW(UbException(
+ UB_EXARGS, "vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+ block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+ UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+
+ if (UbMath::isNaN(phi[REST]) || UbMath::isInfinity(phi[REST]))
+ UB_THROW(UbException(
+ UB_EXARGS, "phi is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+ block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
+ UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+
+ if (UbMath::isNaN(p1) || UbMath::isInfinity(p1))
+ UB_THROW( UbException(UB_EXARGS,"p1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block="+block->toString()+
+ ", node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+
+ data[index++].push_back(phi[REST]);
+ data[index++].push_back(vx1);
+ data[index++].push_back(vx2);
+ data[index++].push_back(vx3);
+ data[index++].push_back(p1);
+ }
+ }
+ }
+ }
+ maxX1 -= 1;
+ maxX2 -= 1;
+ maxX3 -= 1;
+ // cell vector erstellen
+ for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
+ for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
+ for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
+ if ((SWB = nodeNumbers(ix1, ix2, ix3)) >= 0 && (SEB = nodeNumbers(ix1 + 1, ix2, ix3)) >= 0 &&
+ (NEB = nodeNumbers(ix1 + 1, ix2 + 1, ix3)) >= 0 && (NWB = nodeNumbers(ix1, ix2 + 1, ix3)) >= 0 &&
+ (SWT = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0 && (SET = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0 &&
+ (NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0 &&
+ (NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0) {
+ cells.push_back(makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB,
+ (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET,
+ (unsigned int)NET, (unsigned int)NWT));
+ }
+ }
+ }
+ }
+}
+
+LBMReal WriteMultiphaseQuantitiesCoProcessor::gradX1_phi(const LBMReal *const &h)
+{
+ using namespace D3Q27System;
+ LBMReal sum = 0.0;
+ for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ sum += WEIGTH[k] * DX1[k] * h[k];
+ }
+ return 3.0 * sum;
+}
+LBMReal WriteMultiphaseQuantitiesCoProcessor::gradX2_phi(const LBMReal *const &h)
+{
+ using namespace D3Q27System;
+ LBMReal sum = 0.0;
+ for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ sum += WEIGTH[k] * DX2[k] * h[k];
+ }
+ return 3.0 * sum;
+}
+
+LBMReal WriteMultiphaseQuantitiesCoProcessor::gradX3_phi(const LBMReal *const &h)
+{
+ using namespace D3Q27System;
+ LBMReal sum = 0.0;
+ for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ sum += WEIGTH[k] * DX3[k] * h[k];
+ }
+ return 3.0 * sum;
+}
+
+LBMReal WriteMultiphaseQuantitiesCoProcessor::nabla2_phi(const LBMReal *const &h)
+{
+ using namespace D3Q27System;
+ LBMReal sum = 0.0;
+ for (int k = FSTARTDIR; k <= FENDDIR; k++) {
+ sum += WEIGTH[k] * (h[k] - h[REST]);
+ }
+ return 6.0 * sum;
+}
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.h b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4504c1dfccbad377e0bad4bc1aab51989abaff4
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMultiphaseQuantitiesCoProcessor.h
@@ -0,0 +1,104 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids 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.
+//
+// VirtualFluids 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 VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file WriteMultiphaseQuantitiesCoProcessor.h
+//! \ingroup CoProcessors
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#ifndef WriteMultiphaseQuantitiesCoProcessor_H
+#define WriteMultiphaseQuantitiesCoProcessor_H
+
+#include <PointerDefinitions.h>
+#include <string>
+#include <vector>
+
+#include "CoProcessor.h"
+#include "LBMSystem.h"
+#include "UbTuple.h"
+
+class Communicator;
+class Grid3D;
+class UbScheduler;
+class LBMUnitConverter;
+class WbWriter;
+class Block3D;
+
+//! \brief A class writes macroscopic quantities information to a VTK-file
+class WriteMultiphaseQuantitiesCoProcessor : public CoProcessor
+{
+public:
+ WriteMultiphaseQuantitiesCoProcessor();
+ //! \brief Construct WriteMultiphaseQuantitiesCoProcessor object
+ //! \pre The Grid3D and UbScheduler objects must exist
+ //! \param grid is observable Grid3D object
+ //! \param s is UbScheduler object for scheduling of observer
+ //! \param path is path of folder for output
+ //! \param writer is WbWriter object
+ //! \param conv is LBMUnitConverter object
+ //! \param comm is Communicator object
+ WriteMultiphaseQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path,
+ WbWriter *const writer, SPtr<LBMUnitConverter> conv, SPtr<Communicator> comm);
+ ~WriteMultiphaseQuantitiesCoProcessor() override = default;
+
+ void process(double step) override;
+
+protected:
+ //! Collect data for VTK-file
+ //! \param step is a time step
+ void collectData(double step);
+ //! Collect data for VTK-file
+ //! \param block is a time step
+ void addDataMQ(SPtr<Block3D> block);
+ void clearData();
+
+private:
+ void init();
+ std::vector<UbTupleFloat3> nodes;
+ std::vector<UbTupleUInt8> cells;
+ std::vector<std::string> datanames;
+ std::vector<std::vector<double>> data;
+ std::string path;
+ WbWriter *writer;
+ SPtr<LBMUnitConverter> conv;
+ std::vector<std::vector<SPtr<Block3D>>> blockVector;
+ int minInitLevel;
+ int maxInitLevel;
+ int gridRank;
+ SPtr<Communicator> comm;
+
+ LBMReal gradX1_phi(const LBMReal *const &);
+ LBMReal gradX2_phi(const LBMReal *const &);
+ LBMReal gradX3_phi(const LBMReal *const &);
+ LBMReal nabla2_phi(const LBMReal *const &);
+
+ using CalcMacrosFct = void (*)(const LBMReal *const &, LBMReal &, LBMReal &, LBMReal &, LBMReal &);
+ CalcMacrosFct calcMacros;
+};
+
+#endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
index 1bfa01272911c6da84208a9f8484ea6f2be5b78c..f8119ea3557cc528407c560f415c0d7115927a60 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphaseScratchCumulantLBMKernel.cpp
@@ -207,7 +207,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
}
LBMReal collFactorM;
- LBMReal forcingTerm[D3Q27System::ENDF + 1];
+ //LBMReal forcingTerm[D3Q27System::ENDF + 1];
for (int x3 = minX3; x3 < maxX3; x3++) {
for (int x2 = minX2; x2 < maxX2; x2++) {
@@ -286,7 +286,7 @@ void MultiphaseScratchCumulantLBMKernel::calculate(int step)
collFactorM = collFactorL + (collFactorL - collFactorG) * (phi[REST] - phiH) / (phiH - phiL);
- LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
+ //LBMReal mu = 2 * beta * phi[REST] * (phi[REST] - 1) * (2 * phi[REST] - 1) - kappa * nabla2_phi();
//----------- Calculating Macroscopic Values -------------
LBMReal rho = rhoH + rhoToPhi * (phi[REST] - phiH);