diff --git a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
index 1650b9af939195a07996725ac2f1c8b5acd1dc9a..7dff8a88d30317917635c5267c3f75606c9831fb 100644
--- a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
+++ b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
@@ -80,7 +80,7 @@ const real velocity = 9.0;
const real mach = 0.1;
-const uint nodes_per_D = 8;
+const uint nodes_per_D = 32;
std::string path(".");
@@ -88,7 +88,7 @@ std::string simulationName("ActuatorLine");
const uint timeStepOut = 500;
const uint timeStepEnd = 1000;
-
+const float tEnd = 200;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -172,7 +172,7 @@ void multipleLevel(const std::string& configPath)
});
para->setTOut( timeStepOut );
- para->setTEnd( timeStepEnd );
+ para->setTEnd( uint(tEnd/dt) );
para->setIsBodyForce( true );
@@ -199,12 +199,12 @@ void multipleLevel(const std::string& configPath)
ActuatorLine* actuator_line = new ActuatorLine((unsigned int) 3, density, (unsigned int)32, epsilon, turbPos[0], turbPos[1], turbPos[2], D, level, dt, dx);
para->addActuator( actuator_line );
- Probe* probe = new Probe("probe", 100, 100);
+ Probe* probe = new Probe("probe", 100, 500, 100);
std::vector<real> probeCoordsX = {D,2*D,5*D};
std::vector<real> probeCoordsY = {3*D,3*D,3*D};
std::vector<real> probeCoordsZ = {3*D,3*D,3*D};
- // probe->setProbePointsFromList(probeCoordsX,probeCoordsY,probeCoordsZ);
- probe->setProbePointsFromXNormalPlane(2*D, 0.0, 0.0, L_y, L_z, dx, dx);
+ probe->setProbePointsFromList(probeCoordsX,probeCoordsY,probeCoordsZ);
+ // probe->setProbePointsFromXNormalPlane(2*D, 0.0, 0.0, L_y, L_z, dx, dx);
probe->addPostProcessingVariable(PostProcessingVariable::Means);
probe->addPostProcessingVariable(PostProcessingVariable::Variances);
para->addProbe( probe );
diff --git a/src/gpu/VirtualFluids_GPU/Visitor/Probe.cu b/src/gpu/VirtualFluids_GPU/Visitor/Probe.cu
index adf2839d3478da83ffc32e23e91462f9f89df586..f6976c5680e65ec45f6fd5f723c1235e935c8dca 100644
--- a/src/gpu/VirtualFluids_GPU/Visitor/Probe.cu
+++ b/src/gpu/VirtualFluids_GPU/Visitor/Probe.cu
@@ -37,6 +37,51 @@ std::vector<std::string> getPostProcessingVariableNames(PostProcessingVariable v
return varNames;
}
+__device__ void calculateQuantities(uint n, real* quantityArray, bool* quantities, uint* quantityArrayOffsets, uint nPoints, uint node, real vx, real vy, real vz, real rho)
+{
+ //"https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm"
+ // also has extensions for higher order and covariances
+ real inv_n = 1/real(n);
+
+ if(quantities[int(PostProcessingVariable::Means)])
+ {
+ uint arrOff = quantityArrayOffsets[int(PostProcessingVariable::Means)];
+ real vx_m_old = quantityArray[(arrOff+0)*nPoints+node];
+ real vy_m_old = quantityArray[(arrOff+1)*nPoints+node];
+ real vz_m_old = quantityArray[(arrOff+2)*nPoints+node];
+ real rho_m_old = quantityArray[(arrOff+3)*nPoints+node];
+
+ real vx_m_new = ( (n-1)*vx_m_old + vx )*inv_n;
+ real vy_m_new = ( (n-1)*vy_m_old + vy )*inv_n;
+ real vz_m_new = ( (n-1)*vz_m_old + vz )*inv_n;
+ real rho_m_new = ( (n-1)*rho_m_old+ rho )*inv_n;
+
+ quantityArray[(arrOff+0)*nPoints+node] = vx_m_new;
+ quantityArray[(arrOff+1)*nPoints+node] = vy_m_new;
+ quantityArray[(arrOff+2)*nPoints+node] = vz_m_new;
+ quantityArray[(arrOff+3)*nPoints+node] = rho_m_new;
+
+ if(quantities[int(PostProcessingVariable::Variances)])
+ {
+ arrOff = quantityArrayOffsets[int(PostProcessingVariable::Variances)];
+
+ real vx_var_old = quantityArray[(arrOff+0)*nPoints+node];
+ real vy_var_old = quantityArray[(arrOff+1)*nPoints+node];
+ real vz_var_old = quantityArray[(arrOff+2)*nPoints+node];
+ real rho_var_old = quantityArray[(arrOff+3)*nPoints+node];
+
+ real vx_var_new = ( (n-1)*(vx_var_old )+(vx - vx_m_old )*(vx - vx_m_new ) )*inv_n;
+ real vy_var_new = ( (n-1)*(vy_var_old )+(vy - vy_m_old )*(vy - vy_m_new ) )*inv_n;
+ real vz_var_new = ( (n-1)*(vz_var_old )+(vz - vz_m_old )*(vz - vz_m_new ) )*inv_n;
+ real rho_var_new = ( (n-1)*(rho_var_old)+(rho - rho_m_old)*(rho - rho_m_new) )*inv_n;
+
+ quantityArray[(arrOff+0)*nPoints+node] = vx_var_new;
+ quantityArray[(arrOff+1)*nPoints+node] = vy_var_new;
+ quantityArray[(arrOff+2)*nPoints+node] = vz_var_new;
+ quantityArray[(arrOff+3)*nPoints+node] = rho_var_new;
+ }
+ }
+}
__global__ void interpQuantities( uint* pointIndices,
uint nPoints, uint n,
@@ -70,52 +115,13 @@ __global__ void interpQuantities( uint* pointIndices,
real u_interpX, u_interpY, u_interpZ, rho_interp;
- u_interpX = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vx );
- u_interpY = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vy );
- u_interpZ = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vz );
+ u_interpX = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vx );
+ u_interpY = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vy );
+ u_interpZ = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, vz );
rho_interp = trilinearInterpolation( dW, dE, dN, dS, dT, dB, k, ke, kn, kt, kne, kte, ktn, ktne, rho );
- //"https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm"
- // also has extensions for higher order and covariances
- real inv_n = 1/real(n);
- if(quantities[int(PostProcessingVariable::Means)])
- {
- uint arrOff = quantityArrayOffsets[int(PostProcessingVariable::Means)];
- real u_oldX = quantityArray[(arrOff+0)*nPoints+node];
- real u_oldY = quantityArray[(arrOff+1)*nPoints+node];
- real u_oldZ = quantityArray[(arrOff+2)*nPoints+node];
- real rho_old = quantityArray[(arrOff+3)*nPoints+node];
-
- real u_newX = ( (n-1)*u_oldX +u_interpX )*inv_n;
- real u_newY = ( (n-1)*u_oldY +u_interpY )*inv_n;
- real u_newZ = ( (n-1)*u_oldZ +u_interpZ )*inv_n;
- real rho_new = ( (n-1)*rho_old+rho_interp )*inv_n;
-
- quantityArray[(arrOff+0)*nPoints+node] = u_newX;
- quantityArray[(arrOff+1)*nPoints+node] = u_newY;
- quantityArray[(arrOff+2)*nPoints+node] = u_newZ;
- quantityArray[(arrOff+3)*nPoints+node] = rho_new;
-
- if(quantities[int(PostProcessingVariable::Variances)])
- {
- arrOff = quantityArrayOffsets[int(PostProcessingVariable::Variances)];
-
- real u_var_oldX = quantityArray[(arrOff+0)*nPoints+node];
- real u_var_oldY = quantityArray[(arrOff+1)*nPoints+node];
- real u_var_oldZ = quantityArray[(arrOff+2)*nPoints+node];
- real rho_var_old = quantityArray[(arrOff+3)*nPoints+node];
+ calculateQuantities(n, quantityArray, quantities, quantityArrayOffsets, nPoints, node, u_interpX, u_interpY, u_interpZ, rho_interp);
- real u_var_newX = ( (n-1)*(u_var_oldX )+(u_interpX -u_oldX )*(u_interpX -u_newX ) )*inv_n;
- real u_var_newY = ( (n-1)*(u_var_oldY )+(u_interpY -u_oldY )*(u_interpY -u_newY ) )*inv_n;
- real u_var_newZ = ( (n-1)*(u_var_oldZ )+(u_interpZ -u_oldZ )*(u_interpZ -u_newZ ) )*inv_n;
- real rho_var_new = ( (n-1)*(rho_var_old)+(rho_interp-rho_old)*(rho_interp-rho_new) )*inv_n;
-
- quantityArray[(arrOff+0)*nPoints+node] = u_var_newX;
- quantityArray[(arrOff+1)*nPoints+node] = u_var_newY;
- quantityArray[(arrOff+2)*nPoints+node] = u_var_newZ;
- quantityArray[(arrOff+3)*nPoints+node] = rho_var_new;
- }
- }
}
@@ -171,14 +177,6 @@ void Probe::init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager*
std::copy(pointCoordsY_level.begin(), pointCoordsY_level.end(), probeParams[level]->pointCoordsY);
std::copy(pointCoordsZ_level.begin(), pointCoordsZ_level.end(), probeParams[level]->pointCoordsZ);
-
- // probeParams[level]->pointCoordsY = pointCoordsY_level.data();
- // probeParams[level]->pointCoordsZ = pointCoordsZ_level.data();
-
- // for(int i=0; i<probeParams[level]->nPoints; i++)
- // {
- // printf("x %f",probeParams[level]->pointCoordsX[i]);
- // }
// Might have to catch nPoints=0 ?!?!
cudaManager->cudaAllocProbeDistances(this, level);
cudaManager->cudaAllocProbeIndices(this, level);
@@ -225,7 +223,7 @@ void Probe::init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager*
}
-void Probe::visit(Parameter* para, CudaMemoryManager* cudaManager, int level, unsigned int t)
+void Probe::visit(Parameter* para, CudaMemoryManager* cudaManager, int level, uint t)
{
ProbeStruct* probeStruct = this->getProbeStruct(level);
diff --git a/src/gpu/VirtualFluids_GPU/Visitor/Probe.h b/src/gpu/VirtualFluids_GPU/Visitor/Probe.h
index 6b2246fdbd0afd57f7b35435977d8c61839ec15a..54a78d538bbfd418f6d6fc80368e275cf580c5cc 100644
--- a/src/gpu/VirtualFluids_GPU/Visitor/Probe.h
+++ b/src/gpu/VirtualFluids_GPU/Visitor/Probe.h
@@ -46,7 +46,7 @@ public:
}
void init(Parameter* para, GridProvider* gridProvider, CudaMemoryManager* cudaManager);
- void visit(Parameter* para, CudaMemoryManager* cudaManager, int level, unsigned int t);
+ void visit(Parameter* para, CudaMemoryManager* cudaManager, int level, uint t);
void free(Parameter* para, CudaMemoryManager* cudaManager);
ProbeStruct* getProbeStruct(int level){ return this->probeParams[level]; }