From 843dfb84cab19ddbf76a9ee53aec6f8855a0b5b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20Sch=C3=B6nherr?= <schoen@irmb.tu-bs.de> Date: Thu, 9 Nov 2017 13:59:52 +0100 Subject: [PATCH] added output for the turbulent viscosity of the wale model --- .../Calculation/ForceCalculations.cpp | 58 -------- src/VirtualFluids_GPU/GPU/GPU_Interface.h | 1 + src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh | 1 + src/VirtualFluids_GPU/GPU/LBMKernel.cu | 2 + src/VirtualFluids_GPU/GPU/WaleCumulant27.cu | 3 +- .../Interface_OpenFOAM/Interface.cpp | 17 +-- src/VirtualFluids_GPU/LBM/Simulation.cpp | 3 +- .../Output/UnstructuredGridWriter.hpp | 128 ++++++++++++++++++ src/VirtualFluids_GPU/Output/WriteData.cpp | 17 ++- src/VirtualFluids_GPU/Parameter/Parameter.cpp | 25 ++++ src/VirtualFluids_GPU/Parameter/Parameter.h | 5 + 11 files changed, 187 insertions(+), 73 deletions(-) diff --git a/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp b/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp index f670201b5..2bc9b4a17 100644 --- a/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp +++ b/src/VirtualFluids_GPU/Calculation/ForceCalculations.cpp @@ -81,9 +81,6 @@ void ForceCalculations::calcPIDControllerForForce(Parameter* para) tempVeloY /= (double)numberOfElements; tempVeloZ /= (double)numberOfElements; ////////////////////////////////////////////////////////////////// - ////cout debugging - //cout << "tempVeloX = " << tempVeloX << ", tempVeloY = " << tempVeloY << ", tempVeloZ = " << tempVeloZ << endl; - ////////////////////////////////////////////////////////////////// levelVeloAverageX += tempVeloX; levelVeloAverageY += tempVeloY; levelVeloAverageZ += tempVeloZ; @@ -97,10 +94,6 @@ void ForceCalculations::calcPIDControllerForForce(Parameter* para) veloAverageY = levelVeloAverageY / (double)counter; veloAverageZ = levelVeloAverageZ / (double)counter; ////////////////////////////////////////////////////////////////////////// - ////cout debugging - //cout << "veloAverageX = " << veloAverageX << endl; - //cout << "vx1Targed = " << vx1Targed << endl; - ////////////////////////////////////////////////////////////////////////// if (isPID) { //PID-Controller @@ -112,13 +105,6 @@ void ForceCalculations::calcPIDControllerForForce(Parameter* para) y = y / 2.0; } ////////////////////////////////////////////////////////////////////////// - ////cout debugging - //cout << "forceX = " << para->getForcesHost()[0] + y << endl; - //cout << "e = " << e << endl; - //cout << "esum = " << esum << endl; - //cout << "eold = " << eold << endl; - //cout << "y = " << y << endl; - ////////////////////////////////////////////////////////////////////////// para->getForcesDouble()[0] = (para->getForcesDouble()[0] + y); para->getForcesDouble()[1] = (para->getForcesDouble()[1] + y) * 0.0; para->getForcesDouble()[2] = (para->getForcesDouble()[2] + y) * 0.0; @@ -130,47 +116,3 @@ void ForceCalculations::calcPIDControllerForForce(Parameter* para) para->cudaCopyForcingToDevice(); ////////////////////////////////////////////////////////////////////////// } - - - -//ForceCalculations::allocVeloForForcing(Parameter* para) -// { -// ////////////////////////////////////////////////////////////////////////// -// for (int i = 0; i <= para->getFine(); i++) -// { -// int numberOfElements = para->getParH(i)->QSlip.kQ; -// if (numberOfElements > 0) -// { -// para->cudaAllocForceVelo(i,numberOfElements); -// } -// } -// ////////////////////////////////////////////////////////////////////////// -// para->getForcesDouble()[0] = (double)para->getForcesHost()[0]; -// para->getForcesDouble()[1] = (double)para->getForcesHost()[1]; -// para->getForcesDouble()[2] = (double)para->getForcesHost()[2]; -// ////////////////////////////////////////////////////////////////////////// -// } - -//ForceCalculations::printForcing(Parameter* para) -// { -// ////////////////////////////////////////////////////////////////////////// -// //cout << "forcingX = " << para->getForcesHost()[0] << " forcingY = " << para->getForcesHost()[1] << " forcingZ = " << para->getForcesHost()[2] << endl; -// ////////////////////////////////////////////////////////////////////////// -// -// ////////////////////////////////////////////////////////////////////////// -// //set filename -// std::string ffname = para->getFName()+StringUtil::toString<int>(para->getMyID())+"_forcing.txt"; -// const char* fname = ffname.c_str(); -// ////////////////////////////////////////////////////////////////////////// -// //set ofstream -// ofstream ostr; -// ////////////////////////////////////////////////////////////////////////// -// //open file -// ostr.open(fname, std::fstream::app); -// ostr << para->getForcesHost()[0] << " " << para->getForcesHost()[1] << " "<< para->getForcesHost()[2]; -// ostr << endl; -// ////////////////////////////////////////////////////////////////////////// -// //close file -// ostr.close(); -// ////////////////////////////////////////////////////////////////////////// -// } diff --git a/src/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/VirtualFluids_GPU/GPU/GPU_Interface.h index 4542cb7bb..e18b2ca27 100644 --- a/src/VirtualFluids_GPU/GPU/GPU_Interface.h +++ b/src/VirtualFluids_GPU/GPU/GPU_Interface.h @@ -286,6 +286,7 @@ extern "C" void KernelWaleCumOneCompSP27(unsigned int numberOfThreads, doubflo* veloY, doubflo* veloZ, doubflo* DD, + doubflo* turbulentViscosity, int size_Mat, int size_Array, int level, diff --git a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh index 50353d5a6..a2a32bc3f 100644 --- a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh +++ b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh @@ -256,6 +256,7 @@ extern "C" __global__ void LB_Kernel_Wale_Cum_One_Comp_SP_27(doubflo omega, doubflo* veloY, doubflo* veloZ, doubflo* DDStart, + doubflo* turbulentViscosity, int size_Mat, int level, doubflo* forces, diff --git a/src/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/VirtualFluids_GPU/GPU/LBMKernel.cu index 61cba9eb4..32259373b 100644 --- a/src/VirtualFluids_GPU/GPU/LBMKernel.cu +++ b/src/VirtualFluids_GPU/GPU/LBMKernel.cu @@ -913,6 +913,7 @@ extern "C" void KernelWaleCumOneCompSP27(unsigned int numberOfThreads, doubflo* veloY, doubflo* veloZ, doubflo* DD, + doubflo* turbulentViscosity, int size_Mat, int size_Array, int level, @@ -948,6 +949,7 @@ extern "C" void KernelWaleCumOneCompSP27(unsigned int numberOfThreads, veloY, veloZ, DD, + turbulentViscosity, size_Mat, level, forces, diff --git a/src/VirtualFluids_GPU/GPU/WaleCumulant27.cu b/src/VirtualFluids_GPU/GPU/WaleCumulant27.cu index 59f09426f..7730da412 100644 --- a/src/VirtualFluids_GPU/GPU/WaleCumulant27.cu +++ b/src/VirtualFluids_GPU/GPU/WaleCumulant27.cu @@ -14,6 +14,7 @@ extern "C" __global__ void LB_Kernel_Wale_Cum_One_Comp_SP_27(doubflo omega_in, doubflo* veloY, doubflo* veloZ, doubflo* DDStart, + doubflo* turbulentViscosity, int size_Mat, int level, doubflo* forces, @@ -256,7 +257,7 @@ extern "C" __global__ void LB_Kernel_Wale_Cum_One_Comp_SP_27(doubflo omega_in, //nuTurb = rho * powf(delta, two) * powf(SumSd*SumSd, c3o2) / (powf(SumS*SumS, c5o2) + powf(SumSd*SumSd, c5o4)); //nuTurb = powf(delta, two) * powf(SumSd*SumSd, c3o2) / (powf(SumS*SumS, c5o2) + powf(SumSd*SumSd, c5o4)); } - + turbulentViscosity[k] = nuTurb; //////////////////////////////////////////////////////////////////////////////////// //the force be with you { diff --git a/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp b/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp index 009998041..de7812757 100644 --- a/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp +++ b/src/VirtualFluids_GPU/Interface_OpenFOAM/Interface.cpp @@ -40,7 +40,6 @@ void Interface::allocArrays_CoordNeighborGeo(Parameter* para) { //Particles or Wale if (para->getCalcParticle() || para->getUseWale()) { - //cout << "alloc neighborWSB" << endl; neighWSB = new CoordNeighborGeoV(para->getneighborWSB(), binaer, false); } @@ -87,8 +86,9 @@ void Interface::allocArrays_CoordNeighborGeo(Parameter* para) { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// para->cudaAllocCoord(i); para->cudaAllocSP(i); - if (para->getCalcMedian()) para->cudaAllocMedianSP(i); + if (para->getCalcMedian()) para->cudaAllocMedianSP(i); if (para->getCalcParticle() || para->getUseWale()) para->cudaAllocNeighborWSB(i); + if (para->getUseWale()) para->cudaAllocTurbulentViscosity(i); //cout <<"Test 2 " <<endl; //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// coordX->initArrayCoord(para->getParH(i)->coordX_SP, i); @@ -103,12 +103,6 @@ void Interface::allocArrays_CoordNeighborGeo(Parameter* para) { if (para->getCalcParticle() || para->getUseWale()) { neighWSB->initArray(para->getParH(i)->neighborWSB_SP, i); - ////Test - //for (int jj = 1; jj < 100 /*temp*/; jj++) - //{ - // cout << "indexNeighborWSB: " << para->getParH(i)->neighborX_SP[jj] << endl; - // cout << "indexNeighborWSB: " << para->getParH(i)->neighborWSB_SP[jj] << endl; - //} } //cout <<"Test 4 " <<endl; //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -167,6 +161,9 @@ void Interface::allocArrays_CoordNeighborGeo(Parameter* para) { para->getParH(i)->rho_SP_Med[j] = 0.0f; para->getParH(i)->press_SP_Med[j] = 0.0f; } + if (para->getUseWale()) { + para->getParH(i)->turbViscosity[j] = 0.0f; + } } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////Acoustic Test @@ -280,6 +277,10 @@ void Interface::allocArrays_CoordNeighborGeo(Parameter* para) { para->cudaCopyCoord(i); if (para->getCalcMedian()) para->cudaCopyMedianSP(i); if (para->getCalcParticle()) para->cudaCopyNeighborWSB(i); + if (para->getUseWale()) { + para->cudaCopyNeighborWSB(i); + para->cudaCopyTurbulentViscosityHD(i); + } //cout <<"Test 6 " <<endl; //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// } diff --git a/src/VirtualFluids_GPU/LBM/Simulation.cpp b/src/VirtualFluids_GPU/LBM/Simulation.cpp index 429d7248f..18d83ad04 100644 --- a/src/VirtualFluids_GPU/LBM/Simulation.cpp +++ b/src/VirtualFluids_GPU/LBM/Simulation.cpp @@ -558,7 +558,8 @@ void Simulation::run() para->getParD(0)->vx_SP, para->getParD(0)->vy_SP, para->getParD(0)->vz_SP, - para->getParD(0)->d0SP.f[0], + para->getParD(0)->d0SP.f[0], + para->getParD(0)->turbViscosity, para->getParD(0)->size_Mat_SP, para->getParD(0)->size_Array_SP, 0, diff --git a/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp b/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp index a08c18e10..4284a1792 100644 --- a/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp +++ b/src/VirtualFluids_GPU/Output/UnstructuredGridWriter.hpp @@ -314,6 +314,134 @@ namespace UnstrucuredGridWriter + ////////////////////////////////////////////////////////////////////////// + void writeUnstrucuredGridLTwithTurbulentViscosity(Parameter* para, int level, vector<string >& fname) + { + vector< UbTupleFloat3 > nodes; + vector< UbTupleUInt8 > cells; + //vector< UbTupleUInt8 > cells2; + vector< string > nodedatanames; + nodedatanames.push_back("press"); + nodedatanames.push_back("rho"); + nodedatanames.push_back("vx1"); + nodedatanames.push_back("vx2"); + nodedatanames.push_back("vx3"); + nodedatanames.push_back("geo"); + nodedatanames.push_back("turbVis"); + unsigned int number1, number2, number3, number4, number5, number6, number7, number8; + unsigned int dn1, dn2, dn3, dn4, dn5, dn6, dn7, dn8; + bool neighborsFluid; + double vxmax = 0; + unsigned int startpos = 0; + unsigned int endpos = 0; + unsigned int sizeOfNodes = 0; + vector< vector< double > > nodedata(nodedatanames.size()); + + + //printf("\n test for if... \n"); + for (unsigned int part = 0; part < fname.size(); part++) + { + vxmax = 0; + //printf("\n test in if I... \n"); + ////////////////////////////////////////////////////////////////////////// + if (((part + 1)*para->getlimitOfNodesForVTK()) > para->getParH(level)->size_Mat_SP) + { + sizeOfNodes = para->getParH(level)->size_Mat_SP - (part * para->getlimitOfNodesForVTK()); + } + else + { + sizeOfNodes = para->getlimitOfNodesForVTK(); + } + ////////////////////////////////////////////////////////////////////////// + startpos = part * para->getlimitOfNodesForVTK(); + endpos = startpos + sizeOfNodes; + ////////////////////////////////////////////////////////////////////////// + cells.clear(); + nodes.resize(sizeOfNodes); + nodedata[0].resize(sizeOfNodes); + nodedata[1].resize(sizeOfNodes); + nodedata[2].resize(sizeOfNodes); + nodedata[3].resize(sizeOfNodes); + nodedata[4].resize(sizeOfNodes); + nodedata[5].resize(sizeOfNodes); + nodedata[6].resize(sizeOfNodes); + ////////////////////////////////////////////////////////////////////////// + //int counter = 0; + ////////////////////////////////////////////////////////////////////////// + //printf("\n test in if II... \n"); + + for (unsigned int pos = startpos; pos < endpos; pos++) + { + if (/*para->getParH(level)->geoSP[pos] >= GEO_FLUID*/true) + { + ////////////////////////////////////////////////////////////////////////// + double x1 = para->getParH(level)->coordX_SP[pos]; + double x2 = para->getParH(level)->coordY_SP[pos]; + double x3 = para->getParH(level)->coordZ_SP[pos]; + ////////////////////////////////////////////////////////////////////////// + number1 = pos; + dn1 = pos - startpos; + neighborsFluid = true; + ////////////////////////////////////////////////////////////////////////// + nodes[dn1] = (makeUbTuple((float)(x1), (float)(x2), (float)(x3))); + nodedata[0][dn1] = (double)para->getParH(level)->press_SP[pos] / (double)3.0 * (double)para->getDensityRatio() * (double)para->getVelocityRatio() * (double)para->getVelocityRatio(); + nodedata[1][dn1] = (double)para->getParH(level)->rho_SP[pos] / (double)3.0 * (double)para->getDensityRatio() * (double)para->getVelocityRatio() * (double)para->getVelocityRatio(); + nodedata[2][dn1] = (double)para->getParH(level)->vx_SP[pos] * (double)para->getVelocityRatio(); + nodedata[3][dn1] = (double)para->getParH(level)->vy_SP[pos] * (double)para->getVelocityRatio(); + nodedata[4][dn1] = (double)para->getParH(level)->vz_SP[pos] * (double)para->getVelocityRatio(); + nodedata[5][dn1] = (double)para->getParH(level)->geoSP[pos]; + nodedata[6][dn1] = (double)para->getParH(level)->turbViscosity[pos]; + ////////////////////////////////////////////////////////////////////////// + number2 = para->getParH(level)->neighborX_SP[number1]; + number3 = para->getParH(level)->neighborY_SP[number2]; + number4 = para->getParH(level)->neighborY_SP[number1]; + number5 = para->getParH(level)->neighborZ_SP[number1]; + number6 = para->getParH(level)->neighborZ_SP[number2]; + number7 = para->getParH(level)->neighborZ_SP[number3]; + number8 = para->getParH(level)->neighborZ_SP[number4]; + ////////////////////////////////////////////////////////////////////////// + if (para->getParH(level)->geoSP[number2] < GEO_FLUID || + para->getParH(level)->geoSP[number3] < GEO_FLUID || + para->getParH(level)->geoSP[number4] < GEO_FLUID || + para->getParH(level)->geoSP[number5] < GEO_FLUID || + para->getParH(level)->geoSP[number6] < GEO_FLUID || + para->getParH(level)->geoSP[number7] < GEO_FLUID || + para->getParH(level)->geoSP[number8] < GEO_FLUID) neighborsFluid = false; + ////////////////////////////////////////////////////////////////////////// + if (number2 > endpos || + number3 > endpos || + number4 > endpos || + number5 > endpos || + number6 > endpos || + number7 > endpos || + number8 > endpos) neighborsFluid = false; + ////////////////////////////////////////////////////////////////////////// + dn2 = number2 - startpos; + dn3 = number3 - startpos; + dn4 = number4 - startpos; + dn5 = number5 - startpos; + dn6 = number6 - startpos; + dn7 = number7 - startpos; + dn8 = number8 - startpos; + ////////////////////////////////////////////////////////////////////////// + if (isPeriodicCell(para, level, number2, number1, number3, number5)) + continue; + ////////////////////////////////////////////////////////////////////////// + if (neighborsFluid) cells.push_back(makeUbTuple(dn1, dn2, dn3, dn4, dn5, dn6, dn7, dn8)); + ////////////////////////////////////////////////////////////////////////// + } + } + WbWriterVtkXmlBinary::getInstance()->writeOctsWithNodeData(fname[part], nodes, cells, nodedatanames, nodedata); + //WbWriterVtkXmlBinary::getInstance()->writeNodesWithNodeData(fname[part], nodes, nodedatanames, nodedata); + } + } + ////////////////////////////////////////////////////////////////////////// + + + + + + ////////////////////////////////////////////////////////////////////////// void writeUnstrucuredGridLTConc(Parameter* para, int level, vector<string >& fname) { diff --git a/src/VirtualFluids_GPU/Output/WriteData.cpp b/src/VirtualFluids_GPU/Output/WriteData.cpp index 9bee4b524..3475146af 100644 --- a/src/VirtualFluids_GPU/Output/WriteData.cpp +++ b/src/VirtualFluids_GPU/Output/WriteData.cpp @@ -19,10 +19,8 @@ void writeInit(Parameter* para) //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //copy Data to host para->cudaCopyPrint(lev); - if (para->getCalcMedian()) - { - para->cudaCopyMedianPrint(lev); - } + if (para->getCalcMedian()) para->cudaCopyMedianPrint(lev); + if (para->getUseWale()) para->cudaCopyTurbulentViscosityDH(lev); //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// const unsigned int numberOfParts = para->getParH(lev)->size_Mat_SP / para->getlimitOfNodesForVTK() + 1; std::vector<std::string> fname; @@ -76,7 +74,16 @@ void writeInit(Parameter* para) { - UnstrucuredGridWriter::writeUnstrucuredGridLT(para, lev, fname); + + + if (para->getUseWale()) + { + UnstrucuredGridWriter::writeUnstrucuredGridLTwithTurbulentViscosity(para, lev, fname); + } + else + { + UnstrucuredGridWriter::writeUnstrucuredGridLT(para, lev, fname); + } //Debug //InterfaceDebugWriter::writeInterfaceLinesDebugCF(para); diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/VirtualFluids_GPU/Parameter/Parameter.cpp index 6a9304400..3f024e0e9 100644 --- a/src/VirtualFluids_GPU/Parameter/Parameter.cpp +++ b/src/VirtualFluids_GPU/Parameter/Parameter.cpp @@ -563,6 +563,31 @@ void Parameter::cudaFreeNeighborWSB(int lev) { checkCudaErrors( cudaFreeHost(parH[lev]->neighborWSB_SP)); } +//turbulent viscosity +void Parameter::cudaAllocTurbulentViscosity(int lev) +{ + //Host + checkCudaErrors(cudaMallocHost((void**) &(parH[lev]->turbViscosity), parH[lev]->mem_size_doubflo_SP)); + //Device + checkCudaErrors(cudaMalloc((void**) &(parD[lev]->turbViscosity), parD[lev]->mem_size_doubflo_SP)); + ////////////////////////////////////////////////////////////////////////// + double tmp = (double)parH[lev]->mem_size_int_SP; + setMemsizeGPU(tmp, false); +} +void Parameter::cudaCopyTurbulentViscosityHD(int lev) +{ + //copy host to device + checkCudaErrors(cudaMemcpy(parD[lev]->turbViscosity, parH[lev]->turbViscosity, parH[lev]->mem_size_doubflo_SP, cudaMemcpyHostToDevice)); +} +void Parameter::cudaCopyTurbulentViscosityDH(int lev) +{ + //copy device to host + checkCudaErrors(cudaMemcpy(parH[lev]->turbViscosity, parD[lev]->turbViscosity, parH[lev]->mem_size_doubflo_SP, cudaMemcpyDeviceToHost)); +} +void Parameter::cudaFreeTurbulentViscosity(int lev) +{ + checkCudaErrors(cudaFreeHost(parH[lev]->turbViscosity)); +} //median void Parameter::cudaAllocMedianSP(int lev) { diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.h b/src/VirtualFluids_GPU/Parameter/Parameter.h index a59e009ba..f7360b383 100644 --- a/src/VirtualFluids_GPU/Parameter/Parameter.h +++ b/src/VirtualFluids_GPU/Parameter/Parameter.h @@ -301,6 +301,11 @@ public: void cudaCopyNeighborWSB(int lev); void cudaFreeNeighborWSB(int lev); + void cudaAllocTurbulentViscosity(int lev); + void cudaCopyTurbulentViscosityHD(int lev); + void cudaCopyTurbulentViscosityDH(int lev); + void cudaFreeTurbulentViscosity(int lev); + void cudaAllocMedianSP(int lev); void cudaCopyMedianSP(int lev); void cudaFreeMedianSP(int lev); -- GitLab