Skip to content
Snippets Groups Projects
Commit 56facbc5 authored by LEGOLAS\lenz's avatar LEGOLAS\lenz
Browse files

implements infrastructre for flux BCs

parent 386504ed
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 258 additions and 52 deletions
src/GksGpu/BoundaryConditions/BoundaryCondition.cpp
+ 15
− 3
View file @ 56facbc5
......@@ -54,7 +54,13 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
if( cell.type != STOPPER_OUT_OF_GRID && cell.type != STOPPER_OUT_OF_GRID_BOUNDARY ) continue;
for( uint idx = 0; idx < 27; idx++ )
// look in all directions
uint maximalSearchDirection = 27;
// in case of Flux BC look only at face neighbors
if( this->isFluxBC() ) maximalSearchDirection = 6;
for( uint idx = 0; idx < maximalSearchDirection; idx++ )
{
uint neighborCellIdx = cell.cellToCell[ idx ];
......@@ -65,7 +71,7 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
bool neighborCellIsFluid = neighborCell.type != STOPPER_OUT_OF_GRID &&
neighborCell.type != STOPPER_OUT_OF_GRID_BOUNDARY;
bool neighborCellIsValidGhostCell = allowGhostCells && !boundaryFinder( neighborCell.cellCenter );
bool neighborCellIsValidGhostCell = !this->isFluxBC() && allowGhostCells && !boundaryFinder( neighborCell.cellCenter );
if( neighborCellIsFluid || neighborCellIsValidGhostCell )
{
......@@ -79,7 +85,8 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
secondCells.push_back( neighborCell.cellToCell[ idx ] );
}
cell.isWall = this->isWall();
cell.isWall = this->isWall();
cell.isFluxBC = this->isFluxBC();
break;
}
......@@ -100,6 +107,11 @@ void BoundaryCondition::findBoundaryCells(GksMeshAdapter & adapter, bool allowGh
this->myAllocator->allocateMemory( shared_from_this(), ghostCells, domainCells, secondCells );
}
bool BoundaryCondition::isFluxBC()
{
return false;
}
bool BoundaryCondition::secondCellsNeeded()
{
return false;
......
src/GksGpu/BoundaryConditions/BoundaryCondition.h
+ 2
− 0
View file @ 56facbc5
......@@ -43,6 +43,8 @@ struct VF_PUBLIC BoundaryCondition : virtual public BoundaryConditionStruct, pub
virtual bool isWall() = 0;
virtual bool isFluxBC();
virtual bool secondCellsNeeded();
virtual void runBoundaryConditionKernel( const SPtr<DataBase> dataBase,
......
src/GksGpu/BoundaryConditions/InflowComplete.cu
+ 202
− 36
View file @ 56facbc5
......@@ -2,6 +2,7 @@
#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>
#include <cuda.h>
#include <cuda_runtime.h>
......@@ -20,6 +21,12 @@
#include "FlowStateData/FlowStateDataConversion.cuh"
#include "FlowStateData/AccessDeviceData.cuh"
#include "FluxComputation/Moments.cuh"
#include "FluxComputation/ApplyFlux.cuh"
#include "FluxComputation/Transformation.cuh"
#include "FluxComputation/AssembleFlux.cuh"
#include "FluxComputation/ExpansionCoefficients.cuh"
#include "CudaUtility/CudaRunKernel.hpp"
//////////////////////////////////////////////////////////////////////////
......@@ -85,53 +92,207 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
uint ghostCellIdx = boundaryCondition.ghostCells [ startIndex + index ];
uint domainCellIdx = boundaryCondition.domainCells[ startIndex + index ];
PrimitiveVariables ghostCellPrim;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
if( true )
{
PrimitiveVariables domainCellPrim;
PrimitiveVariables ghostCellPrim;
{
PrimitiveVariables domainCellPrim;
{
ConservedVariables domainCellData;
readCellData(domainCellIdx, dataBase, domainCellData);
domainCellPrim = toPrimitiveVariables(domainCellData, parameters.K);
}
// ghostCellPrim.rho = two * boundaryCondition.prim.rho - domainCellPrim.rho;
// ghostCellPrim.U = two * boundaryCondition.prim.U - domainCellPrim.U;
// ghostCellPrim.V = two * boundaryCondition.prim.V - domainCellPrim.V;
// ghostCellPrim.W = two * boundaryCondition.prim.W - domainCellPrim.W;
ghostCellPrim.lambda = /*two * boundaryCondition.prim.lambda -*/ domainCellPrim.lambda;
//#ifdef USE_PASSIVE_SCALAR
// ghostCellPrim.S_1 = two * boundaryCondition.prim.S_1 - domainCellPrim.S_1;
// ghostCellPrim.S_2 = two * boundaryCondition.prim.S_2 - domainCellPrim.S_2;
//#endif // USE_PASSIVE_SCALAR
ghostCellPrim.rho = boundaryCondition.prim.rho;
ghostCellPrim.U = two * boundaryCondition.prim.U - domainCellPrim.U;
ghostCellPrim.V = two * boundaryCondition.prim.V - domainCellPrim.V;
//ghostCellPrim.W = two * boundaryCondition.prim.W - domainCellPrim.W;
ghostCellPrim.W = boundaryCondition.prim.W;
//ghostCellPrim.lambda = boundaryCondition.prim.lambda;
#ifdef USE_PASSIVE_SCALAR
ghostCellPrim.S_1 = boundaryCondition.prim.S_1;
ghostCellPrim.S_2 = boundaryCondition.prim.S_2;
#endif // USE_PASSIVE_SCALAR
real y = dataBase.cellCenter[VEC_Y(ghostCellIdx, dataBase.numberOfCells)];
real x = dataBase.cellCenter[VEC_X(ghostCellIdx, dataBase.numberOfCells)];
real r = sqrt(y*y + x*x);
ghostCellPrim.W *= (one - four*r*r);
#ifdef USE_PASSIVE_SCALAR
ghostCellPrim.S_1 *= (one - four*r*r);
ghostCellPrim.S_2 = boundaryCondition.prim.S_2 - ghostCellPrim.S_1;
#endif // USE_PASSIVE_SCALAR
}
{
ConservedVariables ghostCons = toConservedVariables(ghostCellPrim, parameters.K);
writeCellData(ghostCellIdx, dataBase, ghostCons);
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
if( false )
{
PrimitiveVariables domainCellPrim;
{
ConservedVariables domainCellData;
readCellData( domainCellIdx, dataBase, domainCellData );
domainCellPrim = toPrimitiveVariables( domainCellData, parameters.K );
readCellData(domainCellIdx, dataBase, domainCellData);
domainCellPrim = toPrimitiveVariables(domainCellData, parameters.K);
}
// ghostCellPrim.rho = two * boundaryCondition.prim.rho - domainCellPrim.rho;
// ghostCellPrim.U = two * boundaryCondition.prim.U - domainCellPrim.U;
// ghostCellPrim.V = two * boundaryCondition.prim.V - domainCellPrim.V;
// ghostCellPrim.W = two * boundaryCondition.prim.W - domainCellPrim.W;
ghostCellPrim.lambda = /*two * boundaryCondition.prim.lambda -*/ domainCellPrim.lambda;
//#ifdef USE_PASSIVE_SCALAR
// ghostCellPrim.S_1 = two * boundaryCondition.prim.S_1 - domainCellPrim.S_1;
// ghostCellPrim.S_2 = two * boundaryCondition.prim.S_2 - domainCellPrim.S_2;
//#endif // USE_PASSIVE_SCALAR
ghostCellPrim.rho = boundaryCondition.prim.rho;
ghostCellPrim.U = two * boundaryCondition.prim.U - domainCellPrim.U;
ghostCellPrim.V = two * boundaryCondition.prim.V - domainCellPrim.V;
ghostCellPrim.W = boundaryCondition.prim.W;
//ghostCellPrim.lambda = boundaryCondition.prim.lambda;
#ifdef USE_PASSIVE_SCALAR
ghostCellPrim.S_1 = boundaryCondition.prim.S_1;
ghostCellPrim.S_2 = boundaryCondition.prim.S_2;
#endif // USE_PASSIVE_SCALAR
real y = dataBase.cellCenter[ VEC_Y(ghostCellIdx, dataBase.numberOfCells) ];
real x = dataBase.cellCenter[ VEC_X(ghostCellIdx, dataBase.numberOfCells) ];
real r = sqrt( y*y + x*x );
ghostCellPrim.W *= (one - four*r*r);
#ifdef USE_PASSIVE_SCALAR
ghostCellPrim.S_1 *= (one - four*r*r);
ghostCellPrim.S_2 = one - ghostCellPrim.S_1;
#endif // USE_PASSIVE_SCALAR
real momentU [NUMBER_OF_MOMENTS];
real momentV [NUMBER_OF_MOMENTS];
real momentW [NUMBER_OF_MOMENTS];
real momentXi[NUMBER_OF_MOMENTS];
PrimitiveVariables facePrim = boundaryCondition.prim;
//facePrim.lambda = domainCellPrim.lambda;
transformGlobalToLocal( facePrim, 'z' );
computeMoments(facePrim, parameters.K, momentU, momentV, momentW, momentXi);
ConservedVariables flux;
flux.rho = momentU[0 + 1];
//flux.rhoU = momentU[1 + 1];
//flux.rhoE = c1o2 * ( momentU[2 + 1]
// + momentU[0 + 1] * momentV [2]
// + momentU[0 + 1] * momentW [2]
// + momentU[0 + 1] * momentXi[2] );
flux.rhoE = momentU[0 + 1] * c1o4 * ( parameters.K + five ) / boundaryCondition.prim.lambda;
//////////////////////////////////////////////////////////////////////////
#ifdef USE_PASSIVE_SCALAR
flux.rhoS_1 = flux.rho * boundaryCondition.prim.S_1;
flux.rhoS_2 = flux.rho * boundaryCondition.prim.S_2;
#endif // USE_PASSIVE_SCALAR
flux = ( parameters.dt * parameters.dx * parameters.dx * facePrim.rho ) * flux;
transformLocalToGlobal( flux, 'z' );
applyFluxToPosCell(dataBase, domainCellIdx, flux, 'z', parameters.dt);
return;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
if( false )
{
ConservedVariables ghostCons = toConservedVariables( ghostCellPrim, parameters.K );
PrimitiveVariables domainCellPrim;
{
ConservedVariables domainCellData;
readCellData(domainCellIdx, dataBase, domainCellData);
domainCellPrim = toPrimitiveVariables(domainCellData, parameters.K);
}
PrimitiveVariables facePrim = boundaryCondition.prim;
//////////////////////////////////////////////////////////////////////////
real ax[LENGTH_CELL_DATA];
real ay[LENGTH_CELL_DATA];
real az[LENGTH_CELL_DATA];
real at[LENGTH_CELL_DATA];
#pragma unroll
for( uint i = 0; i < LENGTH_CELL_DATA; i++ )
{
ax[i] = zero;
ay[i] = zero;
az[i] = zero;
at[i] = zero;
}
{
ConservedVariables gradN, gradT1, gradT2;
writeCellData( ghostCellIdx, dataBase, ghostCons );
transformGlobalToLocal( gradN , 'z' );
transformGlobalToLocal( gradT1, 'z' );
transformGlobalToLocal( gradT2, 'z' );
transformGlobalToLocal( facePrim, 'z' );
computeExpansionCoefficients(facePrim, gradN , parameters.K, ax);
computeExpansionCoefficients(facePrim, gradT1, parameters.K, ay);
computeExpansionCoefficients(facePrim, gradT2, parameters.K, az);
}
//////////////////////////////////////////////////////////////////////////
{
ConservedVariables flux;
{
real momentU [ NUMBER_OF_MOMENTS ];
real momentV [ NUMBER_OF_MOMENTS ];
real momentW [ NUMBER_OF_MOMENTS ];
real momentXi[ NUMBER_OF_MOMENTS ];
computeMoments( facePrim, parameters.K, momentU, momentV, momentW, momentXi );
Vec3 force = parameters.force;
transformGlobalToLocal(force, 'z');
{
ConservedVariables timeGrad;
computeTimeDerivative( facePrim,
momentU,
momentV,
momentW,
momentXi,
ax, ay, az,
force,
timeGrad );
computeExpansionCoefficients( facePrim, timeGrad, parameters.K, at );
}
{
real timeCoefficients[4];
computeTimeCoefficients( facePrim, parameters, timeCoefficients );
real heatFlux;
assembleFlux( facePrim,
momentU, momentV, momentW, momentXi,
ax, ay, az, at,
timeCoefficients,
parameters,
force,
flux,
heatFlux );
transformLocalToGlobal( flux, 'z' );
}
}
applyFluxToPosCell(dataBase, domainCellIdx, flux, 'z', parameters.dt);
applyFluxToNegCell(dataBase, ghostCellIdx , flux, 'z', parameters.dt);
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
}
InflowComplete::InflowComplete(SPtr<DataBase> dataBase, PrimitiveVariables prim)
......@@ -145,6 +306,11 @@ bool InflowComplete::isWall()
return false;
}
bool InflowComplete::isFluxBC()
{
return false;
}
bool InflowComplete::secondCellsNeeded()
{
return false;
......
src/GksGpu/BoundaryConditions/InflowComplete.h
+ 2
− 0
View file @ 56facbc5
......@@ -32,6 +32,8 @@ struct VF_PUBLIC InflowComplete : public BoundaryCondition //, public Isothermal
virtual bool isWall() override;
virtual bool isFluxBC() override;
virtual bool secondCellsNeeded() override;
virtual void runBoundaryConditionKernel(const SPtr<DataBase> dataBase,
......
src/GksGpu/BoundaryConditions/Open.cu
+ 10
− 2
View file @ 56facbc5
......@@ -116,8 +116,8 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
//////////////////////////////////////////////////////////////////////////
//if( sign > zero )
if( p2 > p1 )
if( sign < zero )
//if( p2 > p1 )
ghostCellData = domainCellData;
else
{
......@@ -139,6 +139,14 @@ __host__ __device__ inline void boundaryConditionFunction(const DataBaseStruct&
ghostCellData = toConservedVariables(ghostCellPrim, parameters.K);
}
//////////////////////////////////////////////////////////////////////////
//ghostCellData = two * domainCellData + ( - one ) * secondCellData;
//ghostCellData = domainCellData;
//////////////////////////////////////////////////////////////////////////
writeCellData(ghostCellIdx, dataBase, ghostCellData);
}
......
src/GksGpu/CellProperties/CellProperties.cuh
+ 1
− 1
View file @ 56facbc5
......@@ -14,7 +14,7 @@
#define CELL_PROPERTIES_GHOST (1u)
#define CELL_PROPERTIES_WALL (2u)
#define CELL_PROPERTIES_FINE_GHOST (4u)
#define CELL_PROPERTIES_3 (8u)
#define CELL_PROPERTIES_IS_FLUX_BC (8u)
#define CELL_PROPERTIES_4 (16u)
#define CELL_PROPERTIES_5 (32u)
#define CELL_PROPERTIES_6 (64u)
......
src/GksGpu/DataBase/DataBaseAllocatorCPU.cpp
+ 2
− 0
View file @ 56facbc5
......@@ -105,6 +105,8 @@ void DataBaseAllocatorCPU::copyMesh(SPtr<DataBase> dataBase, GksMeshAdapter & ad
dataBase->cellProperties[ cellIdx ] = CELL_PROPERTIES_DEFAULT;
if( adapter.cells[ cellIdx ].isWall )
setCellProperties( dataBase->cellProperties[ cellIdx ], CELL_PROPERTIES_WALL );
if( adapter.cells[ cellIdx ].isFluxBC )
setCellProperties( dataBase->cellProperties[ cellIdx ], CELL_PROPERTIES_IS_FLUX_BC );
if( adapter.cells[ cellIdx ].isGhostCell )
setCellProperties( dataBase->cellProperties[ cellIdx ], CELL_PROPERTIES_GHOST );
if( adapter.cells[ cellIdx ].isFineGhostCell() )
......
src/GksGpu/DataBase/DataBaseAllocatorGPU.cpp
+ 2
− 0
View file @ 56facbc5
......@@ -126,6 +126,8 @@ void DataBaseAllocatorGPU::copyMesh(SPtr<DataBase> dataBase, GksMeshAdapter & ad
cellPropertiesBuffer[ cellIdx ] = CELL_PROPERTIES_DEFAULT;
if( adapter.cells[ cellIdx ].isWall )
setCellProperties( cellPropertiesBuffer[ cellIdx ], CELL_PROPERTIES_WALL );
if( adapter.cells[ cellIdx ].isFluxBC )
setCellProperties( cellPropertiesBuffer[ cellIdx ], CELL_PROPERTIES_IS_FLUX_BC );
if( adapter.cells[ cellIdx ].isGhostCell )
setCellProperties( cellPropertiesBuffer[ cellIdx ], CELL_PROPERTIES_GHOST );
if( adapter.cells[ cellIdx ].isFineGhostCell() )
......
src/GksGpu/FluxComputation/FluxComputation.cu
+ 4
− 2
View file @ 56facbc5
......@@ -190,8 +190,6 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
{
ConservedVariables flux;
ConservedVariables faceCons;
{
real momentU [ NUMBER_OF_MOMENTS ];
real momentV [ NUMBER_OF_MOMENTS ];
......@@ -244,6 +242,10 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
CellProperties negCellProperties = dataBase.cellProperties[ negCellIdx ];
CellProperties posCellProperties = dataBase.cellProperties[ posCellIdx ];
//if( isCellProperties( negCellProperties, CELL_PROPERTIES_IS_FLUX_BC ) ||
// isCellProperties( posCellProperties, CELL_PROPERTIES_IS_FLUX_BC ) )
// return;
if( isCellProperties( negCellProperties, CELL_PROPERTIES_WALL ) ||
isCellProperties( posCellProperties, CELL_PROPERTIES_WALL ) )
{
......
src/GksMeshAdapter/MeshCell.cpp
+ 2
− 0
View file @ 56facbc5
......@@ -21,6 +21,8 @@ MeshCell::MeshCell(){
isGhostCell = false;
isWall = false;
isFluxBC = false;
}
bool MeshCell::isCoarseGhostCell()
......
src/GksMeshAdapter/MeshCell.h
+ 2
− 0
View file @ 56facbc5
......@@ -60,6 +60,8 @@ struct VF_PUBLIC MeshCell{
bool isWall;
bool isFluxBC;
char type;
MeshCell();
......
targets/apps/GKS/Flame7cm/Flame7cm.cpp
+ 14
− 8
View file @ 56facbc5
......@@ -4,6 +4,7 @@
#include <math.h>
#include <string>
#include <iostream>
#include <iomanip>
#include <exception>
#include <fstream>
#include <memory>
......@@ -42,6 +43,7 @@
#include "GksGpu/BoundaryConditions/PassiveScalarDiriclet.h"
#include "GksGpu/BoundaryConditions/InflowComplete.h"
#include "GksGpu/BoundaryConditions/Open.h"
#include "GksGpu/BoundaryConditions/Inflow.h"
#include "GksGpu/Interface/Interface.h"
#include "GksGpu/TimeStepping/NestedTimeStep.h"
......@@ -187,14 +189,17 @@ void thermalCavity( std::string path, std::string simulationName )
SPtr<BoundaryCondition> bcMX = std::make_shared<Open>( dataBase, prim );
SPtr<BoundaryCondition> bcPX = std::make_shared<Open>( dataBase, prim );
//SPtr<BoundaryCondition> bcMX_2 = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
//SPtr<BoundaryCondition> bcPX_2 = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
//SPtr<BoundaryCondition> bcMX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), false );
//SPtr<BoundaryCondition> bcPX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), false );
SPtr<BoundaryCondition> bcMX_2 = std::make_shared<Inflow>( dataBase, Vec3( 100.0*U, 0.0, 0.0), prim.lambda, rho, 1.0, 0.0, 0.0 );
SPtr<BoundaryCondition> bcPX_2 = std::make_shared<Inflow>( dataBase, Vec3(-100.0*U, 0.0, 0.0), prim.lambda, rho, 1.0, 0.0, 0.0 );
bcMX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L; } );
bcPX->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x > 0.5*L; } );
//bcMX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L && center.z > 0.5; } );
//bcPX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x > 0.5*L && center.z > 0.5; } );
bcMX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x < -0.5*L && center.z < 0.25*H; } );
bcPX_2->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ return center.x > 0.5*L && center.z < 0.25*H; } );
//////////////////////////////////////////////////////////////////////////
......@@ -227,19 +232,20 @@ void thermalCavity( std::string path, std::string simulationName )
SPtr<BoundaryCondition> bcPZ = std::make_shared<Open>( dataBase, prim );
bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0; } );
bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0 /*&& std::sqrt(center.x*center.x + center.y*center.y) >= 0.5*0.071*/; } );
bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z > H ; } );
//////////////////////////////////////////////////////////////////////////
//SPtr<BoundaryCondition> burner = std::make_shared<IsothermalWall>( dataBase, Vec3(0.0, 0.0, 0.0), 0.5*prim.lambda, 0.0, true );
SPtr<BoundaryCondition> burner = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, U, prim.lambda, 1.0, 0.0) );
SPtr<BoundaryCondition> burner = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, U, prim.lambda, 1.0, 1.0) );
//SPtr<BoundaryCondition> burner = std::make_shared<InflowComplete>( dataBase, PrimitiveVariables(rho, 0.0, 0.0, 10.0 * U, prim.lambda, 0.0, 0.0) );
burner->findBoundaryCells( meshAdapter, false, [&](Vec3 center){
return center.z < 0.0 && std::sqrt(center.x*center.x + center.y*center.y) < 0.5*0.071;
//return center.z < 0.0 && std::sqrt(center.x*center.x + center.y*center.y) < 0.5*0.071;
return center.z < 0.0 && std::sqrt(center.x*center.x) < 0.5*0.071;
} );
//////////////////////////////////////////////////////////////////////////
......@@ -331,7 +337,7 @@ void thermalCavity( std::string path, std::string simulationName )
if(
//( iter >= 2000 && iter % 100 == 0 ) ||
( iter % 1000 == 0 )
( iter % 10 == 0 )
)
{
for( uint level = 0; level < dataBase->numberOfLevels; level++ )
......
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