MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
Classes | Public Types | Public Member Functions | Static Public Attributes | Protected Attributes | Friends | List of all members
StructuredBndryCnd2D< isRans > Class Template Reference

Class for the 2D stuctured boundary conditions. More...

#include <fvstructuredbndrycnd2d.h>

Inheritance diagram for StructuredBndryCnd2D< isRans >:
[legend]
Collaboration diagram for StructuredBndryCnd2D< isRans >:
[legend]

Classes

struct  comp
 

Public Types

using Timers = maia::structured::Timers_
 
- Public Types inherited from StructuredBndryCnd< 2 >
typedef void(StructuredBndryCnd::* BndryCndHandler) (MInt)
 

Public Member Functions

 StructuredBndryCnd2D (FvStructuredSolver< 2 > *solver, StructuredGrid< 2 > *grid)
 
 ~StructuredBndryCnd2D ()
 
void correctBndryCndIndices ()
 
void correctWallDistanceAtBoundary (MInt)
 
void bc1000 (MInt)
 
template<RansMethod ransMethod>
void bc1000_ (MInt)
 
void bc1001 (MInt)
 
void bc1003 (MInt)
 
void bc1004 (MInt)
 
void bc2001 (MInt)
 Subsonic Inflow <== tfs2001. More...
 
void bc2003 (MInt)
 Subsonic in/outflow simple. More...
 
void bc2004 (MInt)
 Subsonic Outflow, not really non-reflecting for face 0,1,3. More...
 
void bc2002 (MInt)
 Supersonic Inflow. More...
 
void bc2005 (MInt)
 Supersonic Outflow,. More...
 
void bc2006 (MInt)
 Characteristic in/outflow boundary for zero velocities. More...
 
void bc2007 (MInt)
 Subsonic Outflow extrapolate all but pressure, prescribe p8. More...
 
void bc2021 (MInt)
 
void bc2199 (MInt)
 
void bc2402 (MInt)
 
void bc2510 (MInt)
 
void bc2511 (MInt)
 
void bc2600 (MInt)
 Prescribe given profile BC. More...
 
void bc2999 (MInt)
 Blasius bl inflow boundary condition. More...
 
void bc3000 (MInt)
 Symmetry plane BC. More...
 
virtual void bc6002 (MInt) override
 
void initBc1000 (MInt)
 
void initBc1001 (MInt)
 
void initBc1003 (MInt)
 
void initBc1004 (MInt)
 
void initBc2001 (MInt)
 
void initBc2002 (MInt)
 
void initBc2003 (MInt)
 
void initBc2004 (MInt)
 
void initBc2005 (MInt)
 
void initBc2006 (MInt)
 
void initBc2007 (MInt)
 
void initBc2021 (MInt)
 
void initBc2199 (MInt)
 
void initBc2402 (MInt)
 Channel flow BC. More...
 
void initBc2500 (MInt)
 
void initBc2501 (MInt)
 
void initBc2510 (MInt)
 Rescaling inflow. More...
 
void initBc2511 (MInt)
 
void initBc2600 (MInt)
 Prescribe profile BC. More...
 
void initBc3000 (MInt)
 
void initBc6002 (MInt) override
 
virtual void computeFrictionPressureCoef (MBool computePower) override
 
virtual void computeFrictionCoef () override
 
template<MBool calcCp, MBool calcCf, MBool calcIntegrals>
void computeFrictionPressureCoef_ (const MBool auxDataWindow=false, const MBool computePower=false)
 
template<MBool calcCp, MBool calcCf, MBool interface>
void calc_cp_cf (const MInt, const MInt, const MInt, const MInt, MFloat(&)[calcCp+nDim *calcCf])
 
virtual void distributeWallAndFPProperties () override
 
void distributeMapProperties (const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, const std::vector< MInt > &, const std::vector< MInt > &, const std::map< MInt, std::tuple< MInt, MInt, MFloat > > &cellId2recvCell, const std::vector< MInt > &, MFloat *const)
 
void readAndDistributeSpongeCoordinates ()
 
void updateSpongeLayer ()
 
void computeWallDistances ()
 
virtual void computeLocalWallDistances () override
 
void computeDistance2Map (const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, MFloat *const, std::vector< std::pair< MInt, MInt > > &, std::vector< MInt > &, std::vector< MFloat > &)
 Compute shortest distance to given set of maps. More...
 
template<typename T = comp<MFloat>>
void getCloserMap (const MFloat *const, std::vector< std::pair< MInt, MInt > > &, const MFloat *const, std::vector< std::pair< MInt, MInt > > &, MFloat *const, T comparator={})
 
void setUpNearMapComm (const std::vector< std::pair< MInt, MInt > > &, const std::vector< MInt > &, const std::vector< MFloat > &, std::map< MInt, std::tuple< MInt, MInt, MFloat > > &, std::vector< MInt > &, std::vector< MInt > &)
 
virtual void computeLocalExtendedDistancesAndSetComm () override
 
void modifyFPDistance (const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, MFloat *const, std::vector< std::pair< MInt, MInt > > &, const std::vector< MFloat > &)
 Modifies the fluid-porous distance. More...
 
MFloat shortestDistanceToLineElement (const MFloat(&)[nDim], const MFloat(&)[nDim], const MFloat(&)[nDim], MFloat &, MFloat &)
 
MInt cellIndex (MInt i, MInt j)
 
MInt getPointIdFromCell (MInt i, MInt j)
 
MInt getPointIdFromPoint (MInt origin, MInt incI, MInt incJ)
 
MFloat pressure (MInt)
 
MFloat temperature (MInt)
 
- Public Member Functions inherited from StructuredBndryCnd< 2 >
 StructuredBndryCnd (FvStructuredSolver< nDim > *solver, StructuredGrid< nDim > *grid)
 
virtual ~StructuredBndryCnd ()
 
void applyNonReflectingBC ()
 
void assignBndryCnds ()
 
void applyDirichletNeumannBC ()
 
virtual void correctBndryCndIndices ()
 
void saveAuxData ()
 
virtual void bc0 (MInt)
 
virtual void bc1000 (MInt)
 
virtual void bc1001 (MInt)
 
virtual void bc1003 (MInt)
 
virtual void bc1004 (MInt)
 
virtual void bc1006 (MInt)
 
virtual void bc1007 (MInt)
 
virtual void bc2001 (MInt)
 
virtual void bc2002 (MInt)
 
virtual void bc2003 (MInt)
 
virtual void bc2004 (MInt)
 
virtual void bc2024 (MInt)
 
virtual void bc2005 (MInt)
 
virtual void bc2006 (MInt)
 
virtual void bc2007 (MInt)
 
virtual void bc2009 (MInt)
 
virtual void bc2020 (MInt)
 
virtual void bc2021 (MInt)
 
virtual void bc2097 (MInt)
 
virtual void bc2099 (MInt)
 
virtual void bc2221 (MInt)
 
virtual void bc2222 (MInt)
 
virtual void bc2199 (MInt)
 
virtual void bc2402 (MInt)
 
virtual void bc3000 (MInt)
 
virtual void bc3001 (MInt)
 
virtual void bc6000 (MInt)
 
virtual void bc6002 (MInt)
 
virtual void bc2012 (MInt)
 
virtual void bc2014 (MInt)
 
virtual void bc2015 (MInt)
 
virtual void bc7909 (MInt)
 
virtual void bc2013 (MInt)
 
virtual void bc2500 (MInt)
 
virtual void bc2511 (MInt)
 
virtual void bc2510 (MInt)
 
virtual void bc2501 (MInt)
 
virtual void bc2600 (MInt)
 
virtual void bc2601 (MInt)
 
virtual void bc2700 (MInt)
 
virtual void bc2730 (MInt)
 
virtual void bc2888 (MInt)
 
virtual void bc2999 (MInt)
 
virtual void bc2900 (MInt)
 
virtual void initBc0 (MInt)
 
virtual void initBc1000 (MInt)
 
virtual void initBc1001 (MInt)
 
virtual void initBc1003 (MInt)
 
virtual void initBc1004 (MInt)
 
virtual void initBc1006 (MInt)
 
virtual void initBc1007 (MInt)
 
virtual void initBc2001 (MInt)
 
virtual void initBc2003 (MInt)
 
virtual void initBc2004 (MInt)
 
virtual void initBc2024 (MInt)
 
virtual void initBc2002 (MInt)
 
virtual void initBc2005 (MInt)
 
virtual void initBc2006 (MInt)
 
virtual void initBc2007 (MInt)
 
virtual void initBc2009 (MInt)
 
virtual void initBc2020 (MInt)
 
virtual void initBc2021 (MInt)
 
virtual void initBc2097 (MInt)
 
virtual void initBc2099 (MInt)
 
virtual void initBc2221 (MInt)
 
virtual void initBc2222 (MInt)
 
virtual void initBc2199 (MInt)
 
virtual void initBc2402 (MInt)
 
virtual void initBc3000 (MInt)
 
virtual void initBc3001 (MInt)
 
virtual void initBc6000 (MInt)
 
virtual void initBc6002 (MInt)
 
virtual void initBc2012 (MInt)
 
virtual void initBc7909 (MInt)
 
virtual void initBc2013 (MInt)
 
virtual void initBc2014 (MInt)
 
virtual void initBc2015 (MInt)
 
virtual void initBc2500 (MInt)
 
virtual void initBc2501 (MInt)
 
virtual void initBc2510 (MInt)
 
virtual void initBc2600 (MInt)
 
virtual void initBc2601 (MInt)
 
virtual void initBc2700 (MInt)
 
virtual void initBc2730 (MInt)
 
virtual void initBc2888 (MInt)
 
virtual void initBc2999 (MInt)
 
virtual void initBc2900 (MInt)
 
virtual void bc9999 (MInt)
 
virtual void initBc9999 (MInt)
 
virtual void computeWallDistances ()
 
virtual void computeLocalWallDistances ()
 
virtual void computeLocalExtendedDistancesAndSetComm ()
 
virtual void updateSpongeLayer ()
 
virtual void computeFrictionCoef ()
 
virtual void computeFrictionPressureCoef (MBool)=0
 
virtual void distributeWallAndFPProperties ()
 

Static Public Attributes

static constexpr const MInt nDim = 2
 

Protected Attributes

FvStructuredSolver2D * m_solver
 
MInt m_startCommPeriodic
 
MInt m_endCommPeriodic
 
MInt m_startCommChannel
 
MInt m_endCommChannel
 
MInt m_channelInflowRank
 
MInt m_noPeriodicConnections
 
MFloat m_rescalingBLT
 
MFloat m_isothermalWallTemperature
 
MFloat m_bc2021Gradient
 
- Protected Attributes inherited from StructuredBndryCnd< 2 >
std::vector< MInt > m_wallSendCells
 
std::vector< MInt > m_wallSendcounts
 
std::map< MInt, std::tuple< MInt, MInt, MFloat > > m_wallCellId2recvCell
 
std::vector< MInt > m_FPSendCells
 
std::vector< MInt > m_FPSendcounts
 
std::map< MInt, std::tuple< MInt, MInt, MFloat > > m_FPCellId2recvCell
 
std::vector< MInt > m_FPSendCells_porous
 
std::vector< MInt > m_FPSendcounts_porous
 
std::map< MInt, std::tuple< MInt, MInt, MFloat > > m_FPCellId2recvCell_porous
 
MInt m_noCells
 
MFloat m_channelSurfaceIn
 
MFloat m_channelSurfaceOut
 
MInt * m_channelSurfaceIndexMap
 
MInt * m_plenumSurfaceIndexMap
 
MFloat m_plenumSurface
 
MFloat m_sutherlandPlusOne
 
MFloat m_sutherlandConstant
 
MFloat m_UInfinity
 
MFloat m_VInfinity
 
MFloat m_WInfinity
 
MFloat m_PInfinity
 
MFloat m_TInfinity
 
MFloat m_DthInfinity
 
MFloat m_muInfinity
 
MFloat m_DInfinity
 
MFloat m_VVInfinity [3]
 
MFloat m_rhoUInfinity
 
MFloat m_rhoVInfinity
 
MFloat m_rhoWInfinity
 
MFloat m_rhoEInfinity
 
MFloat m_rhoInfinity
 
MFloat m_rhoVVInfinity [3]
 

Friends

class FvStructuredSolver2D
 

Additional Inherited Members

- Public Attributes inherited from StructuredBndryCnd< 2 >
MPI_Comm m_StructuredComm
 
MInt * m_nCells
 
MInt * m_nPoints
 
StructuredCell * m_cells
 
FvStructuredSolver< nDim > * m_solver
 
StructuredGrid< nDim > * m_grid
 
std::unique_ptr< MConservativeVariables< nDim > > & CV
 
std::unique_ptr< MPrimitiveVariables< nDim > > & PV
 
std::unique_ptr< StructuredFQVariables > & FQ
 
MInt m_solverId
 
MInt m_noBndryCndIds
 
MInt m_noGhostLayers
 
std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > & m_physicalBCMap
 
std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > & m_auxDataMap
 
std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > & m_globalStructuredBndryMaps
 
MInt m_noSpongeDomainInfos
 
MInt * m_spongeBcWindowInfo
 
MInt m_spongeLayerType
 
MFloat * m_spongeLayerThickness
 
MFloat * m_sigmaSponge
 
MFloat * m_betaSponge
 
MFloat m_targetDensityFactor
 
MFloat m_sigma
 
BndryCndHandler * nonReflectingBoundaryCondition
 
BndryCndHandler * bndryCndHandler
 
BndryCndHandler * initBndryCndHandler
 
BndryCndHandler * skinFrictionHandler
 

Detailed Description

template<MBool isRans>
class StructuredBndryCnd2D< isRans >

Definition at line 21 of file fvstructuredbndrycnd2d.h.

Member Typedef Documentation

◆ Timers

template<MBool isRans>
using StructuredBndryCnd2D< isRans >::Timers = maia::structured::Timers_

Definition at line 25 of file fvstructuredbndrycnd2d.h.

Constructor & Destructor Documentation

◆ StructuredBndryCnd2D()

template<MBool isRans>
StructuredBndryCnd2D< isRans >::StructuredBndryCnd2D ( FvStructuredSolver< 2 > *  solver,
StructuredGrid< 2 > *  grid 
)

Definition at line 26 of file fvstructuredbndrycnd2d.cpp.

27 : StructuredBndryCnd<2>(solver, grid) {
28 TRACE();
29
30 const MLong oldAllocatedBytes = allocatedBytes();
31
32 m_solver = static_cast<FvStructuredSolver2D*>(solver);
33
34 // read rescaling bc properties
35 if(m_solver->m_rescalingCommGrRoot >= 0) {
36 m_rescalingBLT = 5.95;
37 m_rescalingBLT = Context::getSolverProperty<MFloat>("rescalingBLT", m_solverId, AT_);
38
39 cout << "Rescaling BLT: " << m_rescalingBLT << endl;
40 }
41
42 // compute sponge coefficients for each cell
43 if(m_solver->m_useSponge) {
44 RECORD_TIMER_START(m_solver->timer(Timers::BuildUpSponge));
45 readAndDistributeSpongeCoordinates();
46 RECORD_TIMER_STOP(m_solver->timer(Timers::BuildUpSponge));
47 }
48
49 printAllocatedMemory(oldAllocatedBytes, "StructuredBndryCnd2D", m_StructuredComm);
50}
allocatedBytes
MLong allocatedBytes()
Return the number of allocated bytes.
Definition: alloc.cpp:121
FvStructuredSolver2D
2D structured solver class
Definition: fvstructuredsolver2d.h:21
FvStructuredSolver::timer
MInt timer(const MInt timerId) const
Definition: fvstructuredsolver.cpp:3071
StructuredBndryCnd2D::m_solver
FvStructuredSolver2D * m_solver
Definition: fvstructuredbndrycnd2d.h:139
StructuredBndryCnd
Base class of the structured boundary conditions.
Definition: fvstructuredbndrycnd.h:27
printAllocatedMemory
void printAllocatedMemory(const MLong oldAllocatedBytes, const MString &solverName, const MPI_Comm comm)
Prints currently allocated memory.
Definition: genericobject.cpp:26
MLong
int64_t MLong
Definition: maiatypes.h:64
maia::structured::Timers_::BuildUpSponge
@ BuildUpSponge
Definition: fvstructuredtimers.h:29

◆ ~StructuredBndryCnd2D()

template<MBool isRans>
StructuredBndryCnd2D< isRans >::~StructuredBndryCnd2D

Definition at line 53 of file fvstructuredbndrycnd2d.cpp.

53{}

Member Function Documentation

◆ bc1000()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc1000 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2954 of file fvstructuredbndrycnd2d.cpp.

2954 {
2955 if(isRans) {
2956 const RansMethod ransMethod = m_solver->m_ransMethod;
2957 if(ransMethod == RANS_SA || ransMethod == RANS_SA_DV || ransMethod == RANS_FS) {
2958 bc1000_<RANS_SA>(bcId);
2959 } else if(ransMethod == RANS_KEPSILON) {
2960 bc1000_<RANS_KEPSILON>(bcId);
2961 }
2962 } else
2963 bc1000_<NORANS>(bcId);
2964}
RansMethod
RansMethod
Definition: enums.h:51
RANS_SA_DV
@ RANS_SA_DV
Definition: enums.h:54
RANS_SA
@ RANS_SA
Definition: enums.h:53
RANS_FS
@ RANS_FS
Definition: enums.h:55
RANS_KEPSILON
@ RANS_KEPSILON
Definition: enums.h:58

◆ bc1000_()

template<MBool isRans>
template<RansMethod ransMethod>
void StructuredBndryCnd2D< isRans >::bc1000_ ( MInt  bcId)

Definition at line 2888 of file fvstructuredbndrycnd2d.cpp.

2888 {
2889 TRACE();
2890
2891 const MInt IJ[nDim] = {1, m_nCells[1]};
2892 MInt* start = m_physicalBCMap[bcId]->start1;
2893 MInt* end = m_physicalBCMap[bcId]->end1;
2894
2895 // Here we find out the normal direction of the
2896 // boundary and the tangential direction.
2897 // This way we can make a general formulation of
2898 // the boundary condition
2899 const MInt face = m_physicalBCMap[bcId]->face;
2900 const MInt normalDir = face / 2;
2901 const MInt firstTangentialDir = (normalDir + 1) % nDim;
2902 const MInt normalDirStart = start[normalDir];
2903 const MInt firstTangentialStart = start[firstTangentialDir];
2904 const MInt firstTangentialEnd = end[firstTangentialDir];
2905 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
2906 // determine indices for direction help
2907 const MInt n = (face % 2) * 2 - 1; //-1,+1
2908 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
2909 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
2910 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
2911 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
2912
2913 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
2914 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
2915 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
2916 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
2917 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
2918
2919 const MFloat p1 = m_cells->pvariables[PV->P][cellIdA1];
2920 const MFloat rho1 = m_cells->pvariables[PV->RHO][cellIdA1];
2921 const MFloat u1 = m_cells->pvariables[PV->U][cellIdA1];
2922 const MFloat v1 = m_cells->pvariables[PV->V][cellIdA1];
2923
2924 const MFloat p2 = m_cells->pvariables[PV->P][cellIdA2];
2925 const MFloat rho2 = m_cells->pvariables[PV->RHO][cellIdA2];
2926 const MFloat u2 = m_cells->pvariables[PV->U][cellIdA2];
2927 const MFloat v2 = m_cells->pvariables[PV->V][cellIdA2];
2928
2929 m_cells->pvariables[PV->RHO][cellIdG1] = rho1;
2930 m_cells->pvariables[PV->U][cellIdG1] = -u1;
2931 m_cells->pvariables[PV->V][cellIdG1] = -v1;
2932 m_cells->pvariables[PV->P][cellIdG1] = p1;
2933
2934 m_cells->pvariables[PV->RHO][cellIdG2] = rho2;
2935 m_cells->pvariables[PV->U][cellIdG2] = -u2;
2936 m_cells->pvariables[PV->V][cellIdG2] = -v2;
2937 m_cells->pvariables[PV->P][cellIdG2] = p2;
2938
2939 if(ransMethod == RANS_SA || ransMethod == RANS_KEPSILON) {
2940 // Note: not all k-epsilon models have the same wall BC
2941 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
2942 // For e.g. SA-model the only rans variable is ransVar=nuTilde
2943 const MFloat ransVar1 = m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdA1];
2944 const MFloat ransVar2 = m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdA2];
2945
2946 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdG1] = -ransVar1;
2947 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdG2] = -ransVar2;
2948 }
2949 }
2950 }
2951}
StructuredBndryCnd2D::nDim
static constexpr const MInt nDim
Definition: fvstructuredbndrycnd2d.h:31
StructuredBndryCnd< 2 >::PV
std::unique_ptr< MPrimitiveVariables< nDim > > & PV
Definition: fvstructuredbndrycnd.h:54
StructuredBndryCnd< 2 >::m_physicalBCMap
std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > & m_physicalBCMap
Definition: fvstructuredbndrycnd.h:62
StructuredCell::pvariables
MFloat ** pvariables
Definition: fvstructuredcell.h:18
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52

◆ bc1001()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc1001 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2969 of file fvstructuredbndrycnd2d.cpp.

2969 {
2970 TRACE();
2971
2972 MInt* start = m_physicalBCMap[bcId]->start1;
2973 MInt* end = m_physicalBCMap[bcId]->end1;
2974 switch(m_physicalBCMap[bcId]->face) {
2975 case 2:
2976 case 3: {
2977 MInt cellShift = 0;
2978 if(m_physicalBCMap[bcId]->face == 2) {
2979 cellShift = 2 * m_noGhostLayers - 1;
2980 } else if(m_physicalBCMap[bcId]->face == 3) {
2981 cellShift = 2 * (m_nCells[0] - 1) - 2 * m_noGhostLayers + 1;
2982 }
2983
2984 for(MInt j = start[1]; j < end[1]; j++) {
2985 for(MInt i = start[0]; i < end[0]; i++) {
2986 MInt cellIdA1 = -1, pIJ = -1;
2987 if(m_physicalBCMap[bcId]->face == 2) {
2988 pIJ = getPointIdFromCell(i, m_noGhostLayers);
2989 cellIdA1 = cellIndex(i, m_noGhostLayers);
2990 } else {
2991 pIJ = getPointIdFromCell(i, m_solver->m_nPoints[0] - 3);
2992 cellIdA1 = cellIndex(i, m_nCells[0] - 3);
2993 }
2994 const MFloat x1 = m_grid->m_coordinates[0][pIJ];
2995 const MFloat y1 = m_grid->m_coordinates[1][pIJ];
2996 const MFloat x2 = m_grid->m_coordinates[0][pIJ + 1];
2997 const MFloat y2 = m_grid->m_coordinates[1][pIJ + 1];
2998
2999 const MFloat dydx = (y2 - y1) / (x2 - x1);
3000 const MFloat alpha = -atan(dydx);
3001
3002 const MInt cellId = cellIndex(i, j); // ghost
3003 const MInt cellIdadj = cellIndex(i, cellShift - j); // field
3004
3005 const MFloat rho =
3006 m_cells->pvariables[PV->RHO][cellIdA1]; // apply rho from first active cell to both ghost-cells
3007 const MFloat u = m_cells->pvariables[PV->U][cellIdadj];
3008 const MFloat v = m_cells->pvariables[PV->V][cellIdadj];
3009
3010 const MFloat uPrime = u * cos(alpha) - v * sin(alpha);
3011 const MFloat vPrime = -(u * sin(alpha) + v * cos(alpha));
3012
3013 const MFloat uGC = uPrime * cos(-alpha) - vPrime * sin(-alpha);
3014 const MFloat vGC = uPrime * sin(-alpha) + vPrime * cos(-alpha);
3015
3016 m_cells->pvariables[PV->RHO][cellId] = rho;
3017 m_cells->pvariables[PV->U][cellId] = uGC;
3018 m_cells->pvariables[PV->V][cellId] = vGC;
3019
3020 // apply pressure from first active cell to all cells
3021 m_cells->pvariables[PV->P][cellId] = m_cells->pvariables[PV->P][cellIdA1];
3022
3023 if(isRans) {
3024 m_cells->pvariables[PV->RANS_VAR[0]][cellId] = -m_cells->pvariables[PV->RANS_VAR[0]][cellIdadj];
3025 }
3026 }
3027 }
3028 break;
3029 }
3030 default: {
3031 mTerm(1, AT_, "Face direction not implemented)");
3032 }
3033 }
3034}
StructuredBndryCnd2D::getPointIdFromCell
MInt getPointIdFromCell(MInt i, MInt j)
Definition: fvstructuredbndrycnd2d.cpp:2505
StructuredBndryCnd2D::cellIndex
MInt cellIndex(MInt i, MInt j)
Definition: fvstructuredbndrycnd2d.cpp:2500
StructuredBndryCnd< 2 >::m_grid
StructuredGrid< nDim > * m_grid
Definition: fvstructuredbndrycnd.h:52
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
maia::collector_memory::cellId
void const MInt cellId
Definition: collector.h:239
maia::filter::slope::cos
T cos(const T a, const T b, const T x)
Cosine slope filter.
Definition: filter.h:125

◆ bc1003()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc1003 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3038 of file fvstructuredbndrycnd2d.cpp.

3038 {
3039 const MInt IJ[nDim] = {1, m_nCells[1]};
3040 MInt* start = m_physicalBCMap[bcId]->start1;
3041 MInt* end = m_physicalBCMap[bcId]->end1;
3042
3043 MFloat temp = m_isothermalWallTemperature * PV->TInfinity;
3044 const MFloat gamma = m_solver->m_gamma;
3045
3046 // Here we find out the normal direction of the
3047 // boundary and the tangential direction.
3048 // This way we can make a general formulation of
3049 // the boundary condition
3050 const MInt face = m_physicalBCMap[bcId]->face;
3051 const MInt normalDir = face / 2;
3052 const MInt firstTangentialDir = (normalDir + 1) % nDim;
3053 const MInt normalDirStart = start[normalDir];
3054 const MInt firstTangentialStart = start[firstTangentialDir];
3055 const MInt firstTangentialEnd = end[firstTangentialDir];
3056 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
3057 // determine indices for direction help
3058 const MInt n = (face % 2) * 2 - 1; //-1,+1
3059 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
3060 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
3061 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
3062 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
3063
3064 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
3065 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
3066 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
3067 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
3068 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
3069
3070 const MFloat p1 = m_cells->pvariables[PV->P][cellIdA1];
3071 const MFloat rho1 = m_cells->pvariables[PV->RHO][cellIdA1];
3072 const MFloat u1 = m_cells->pvariables[PV->U][cellIdA1];
3073 const MFloat v1 = m_cells->pvariables[PV->V][cellIdA1];
3074
3075 const MFloat p2 = m_cells->pvariables[PV->P][cellIdA2];
3076 const MFloat rho2 = m_cells->pvariables[PV->RHO][cellIdA2];
3077 const MFloat u2 = m_cells->pvariables[PV->U][cellIdA2];
3078 const MFloat v2 = m_cells->pvariables[PV->V][cellIdA2];
3079
3080 const MFloat rhoWall = p1 * gamma / temp;
3081 const MFloat rhoG1 = F2 * rhoWall - rho1;
3082 const MFloat rhoG2 = F2 * rhoWall - rho2;
3083
3084 m_cells->pvariables[PV->RHO][cellIdG1] = rhoG1;
3085 m_cells->pvariables[PV->U][cellIdG1] = -u1;
3086 m_cells->pvariables[PV->V][cellIdG1] = -v1;
3087 m_cells->pvariables[PV->P][cellIdG1] = p1;
3088
3089 m_cells->pvariables[PV->RHO][cellIdG2] = rhoG2;
3090 m_cells->pvariables[PV->U][cellIdG2] = -u2;
3091 m_cells->pvariables[PV->V][cellIdG2] = -v2;
3092 m_cells->pvariables[PV->P][cellIdG2] = p2;
3093
3094 if(isRans) {
3095 // Note: not all k-epsilon models have the same wall BC
3096 // For e.g. SA-model the only rans variable is ransVar=nuTilde
3097 const MFloat ransVar1 = m_cells->pvariables[PV->RANS_VAR[0]][cellIdA1];
3098 const MFloat ransVar2 = m_cells->pvariables[PV->RANS_VAR[0]][cellIdA2];
3099
3100 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG1] = -ransVar1;
3101 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG2] = -ransVar2;
3102 }
3103 }
3104}

◆ bc1004()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc1004 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3108 of file fvstructuredbndrycnd2d.cpp.

3108 {
3109 TRACE();
3110
3111 MInt* start = m_physicalBCMap[bcId]->start1;
3112 MInt* end = m_physicalBCMap[bcId]->end1;
3113
3114 const MInt IJ[nDim] = {1, m_nCells[1]};
3115 const MInt IJP[nDim] = {1, m_nPoints[1]};
3116
3117 const MInt pp[2][4] = {{0, 0, 0, 1}, {0, 0, 1, 0}};
3118
3119 // Here we find out the normal direction of the
3120 // boundary and the tangential direction.
3121 // This way we can make a general formulation of
3122 // the boundary condition
3123 const MInt face = m_physicalBCMap[bcId]->face;
3124 const MInt normalDir = face / 2;
3125 const MInt firstTangentialDir = (normalDir + 1) % nDim;
3126 const MInt normalDirStart = start[normalDir];
3127 const MInt firstTangentialStart = start[firstTangentialDir];
3128 const MInt firstTangentialEnd = end[firstTangentialDir];
3129 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
3130 const MInt incp[nDim] = {IJP[normalDir], IJP[firstTangentialDir]};
3131 // determine indices for direction help
3132 const MInt n = (face % 2) * 2 - 1; //-1,+1
3133 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
3134 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
3135 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
3136 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
3137
3138 const MInt g1p = normalDirStart + 2 * ((MInt)(0.5 - (0.5 * (MFloat)n))); //+1,0
3139
3140 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
3141 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
3142 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
3143 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
3144 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
3145
3146 const MInt ij = g1p * incp[0] + t1 * incp[1];
3147
3148 const MInt pp1 = getPointIdFromPoint(ij, pp[normalDir][0], pp[normalDir][1]);
3149 const MInt pp2 = getPointIdFromPoint(ij, pp[normalDir][2], pp[normalDir][3]);
3150
3151 // compute the velocity of the surface centroid
3152 MFloat gridVel[2] = {F0, F0};
3153 // MFloat gridAcc[2] = {F0, F0};
3154 // MFloat firstVec[2] = {F0, F0};
3155 // MFloat secondVec[2] = {F0, F0};
3156 // MFloat normalVec[2] = {F0, F0};
3157 for(MInt dim = 0; dim < nDim; dim++) {
3158 // firstVec[dim] = m_grid->m_coordinates[dim][pp2] - m_grid->m_coordinates[dim][pp1];
3159 // secondVec[dim] = m_grid->m_coordinates[dim][pp3] - m_grid->m_coordinates[dim][pp1];
3160 gridVel[dim] = F1B2 * (m_grid->m_velocity[dim][pp1] + m_grid->m_velocity[dim][pp2]);
3161 // gridAcc[dim] = F1B2 * (m_grid->m_acceleration[dim][pp1] + m_grid->m_acceleration[dim][pp2]);
3162 }
3163
3164 const MFloat p1 = m_cells->pvariables[PV->P][cellIdA1];
3165 const MFloat rho1 = m_cells->pvariables[PV->RHO][cellIdA1];
3166 const MFloat u1 = m_cells->pvariables[PV->U][cellIdA1];
3167 const MFloat v1 = m_cells->pvariables[PV->V][cellIdA1];
3168
3169 const MFloat p2 = m_cells->pvariables[PV->P][cellIdA2];
3170 const MFloat rho2 = m_cells->pvariables[PV->RHO][cellIdA2];
3171 const MFloat u2 = m_cells->pvariables[PV->U][cellIdA2];
3172 const MFloat v2 = m_cells->pvariables[PV->V][cellIdA2];
3173
3174 const MFloat pG1 = p1;
3175 const MFloat pG2 = p2;
3176 const MFloat rhoG1 = rho1;
3177 const MFloat rhoG2 = rho2;
3178
3179 const MFloat uG1 = F2 * gridVel[0] - u1;
3180 const MFloat vG1 = F2 * gridVel[1] - v1;
3181
3182 const MFloat uG2 = F2 * gridVel[0] - u2;
3183 const MFloat vG2 = F2 * gridVel[1] - v2;
3184
3185 m_cells->pvariables[PV->RHO][cellIdG1] = rhoG1;
3186 m_cells->pvariables[PV->U][cellIdG1] = uG1;
3187 m_cells->pvariables[PV->V][cellIdG1] = vG1;
3188 m_cells->pvariables[PV->P][cellIdG1] = pG1;
3189
3190 m_cells->pvariables[PV->RHO][cellIdG2] = rhoG2;
3191 m_cells->pvariables[PV->U][cellIdG2] = uG2;
3192 m_cells->pvariables[PV->V][cellIdG2] = vG2;
3193 m_cells->pvariables[PV->P][cellIdG2] = pG2;
3194 }
3195}
StructuredBndryCnd2D::getPointIdFromPoint
MInt getPointIdFromPoint(MInt origin, MInt incI, MInt incJ)
Definition: fvstructuredbndrycnd2d.cpp:2510

◆ bc2001()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2001 ( MInt  bcId)
virtual

rho=rho_inf, u=u_inf, v=v_inf, w=w_inf, dp/dn=0

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3203 of file fvstructuredbndrycnd2d.cpp.

3203 {
3204 TRACE();
3205
3206 // implemented for i-direction only for the moment
3207 MInt* start = m_physicalBCMap[bcId]->start1;
3208 MInt* end = m_physicalBCMap[bcId]->end1;
3209 switch(m_physicalBCMap[bcId]->face) {
3210 case 0: {
3211 MInt cellId = -1;
3212 MInt cellIdadj = -1;
3213 for(MInt j = start[1]; j < end[1]; j++) {
3214 for(MInt i = start[0]; i < end[0]; i++) {
3215 cellId = cellIndex(m_noGhostLayers - 1 - i, j);
3216 cellIdadj = cellIndex(m_noGhostLayers - i, j);
3217
3218 m_cells->pvariables[PV->RHO][cellId] = CV->rhoInfinity;
3219 m_cells->pvariables[PV->U][cellId] = PV->UInfinity;
3220 m_cells->pvariables[PV->V][cellId] = PV->VInfinity;
3221 m_cells->pvariables[PV->P][cellId] = m_cells->pvariables[PV->P][cellIdadj];
3222
3223 if(isRans) {
3224 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3225 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] = PV->ransInfinity[ransVarId];
3226 }
3227 }
3228 }
3229 }
3230 break;
3231 }
3232 case 2: {
3233 MInt cellId = -1;
3234 MInt cellIdadj = -1;
3235 for(MInt j = start[1]; j < end[1]; j++) {
3236 for(MInt i = start[0]; i < end[0]; i++) {
3237 cellId = cellIndex(i, m_noGhostLayers - j - 1); // ghost
3238 cellIdadj = cellIndex(i, m_noGhostLayers - j); // field
3239
3240 m_cells->pvariables[PV->RHO][cellId] = CV->rhoInfinity;
3241 m_cells->pvariables[PV->U][cellId] = PV->UInfinity;
3242 m_cells->pvariables[PV->V][cellId] = PV->VInfinity;
3243 m_cells->pvariables[PV->P][cellId] = m_cells->pvariables[PV->P][cellIdadj];
3244
3245 if(isRans) {
3246 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3247 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] = PV->ransInfinity[ransVarId];
3248 }
3249 }
3250 }
3251 }
3252 break;
3253 }
3254 default: {
3255 mTerm(1, AT_, "Face direction not implemented)");
3256 }
3257 }
3258}
StructuredBndryCnd< 2 >::CV
std::unique_ptr< MConservativeVariables< nDim > > & CV
Definition: fvstructuredbndrycnd.h:53

◆ bc2002()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2002 ( MInt  bcId)
virtual

rho=rho, u=u, v=v, w=w, p=p

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3367 of file fvstructuredbndrycnd2d.cpp.

3367 {
3368 MInt* start = m_physicalBCMap[bcId]->start1;
3369 MInt* end = m_physicalBCMap[bcId]->end1;
3370 switch(m_physicalBCMap[bcId]->face) {
3371 case 1: {
3372 for(MInt j = start[1]; j < end[1]; j++) {
3373 for(MInt i = start[0]; i < end[0]; i++) {
3374 MInt cellId = cellIndex(i, j);
3375 m_cells->pvariables[PV->RHO][cellId] = CV->rhoInfinity;
3376 m_cells->pvariables[PV->U][cellId] = PV->UInfinity;
3377 m_cells->pvariables[PV->V][cellId] = PV->VInfinity;
3378 m_cells->pvariables[PV->P][cellId] = PV->PInfinity;
3379 if(isRans) {
3380 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3381 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] = PV->ransInfinity[ransVarId];
3382 }
3383 }
3384 }
3385 }
3386 break;
3387 }
3388 default: {
3389 // do nothing
3390 }
3391 }
3392}

◆ bc2003()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2003 ( MInt  bcId)
virtual

Simple extrapolation/prescription depending on in/outflow condition

Author
Marian Albers

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3403 of file fvstructuredbndrycnd2d.cpp.

3403 {
3404 const MInt IJ[nDim] = {1, m_nCells[1]};
3405 MInt* start = m_physicalBCMap[bcId]->start1;
3406 MInt* end = m_physicalBCMap[bcId]->end1;
3407
3408 // Here we find out the normal direction of the
3409 // boundary and the tangential direction.
3410 // This way we can make a general formulation of
3411 // the boundary condition
3412 const MInt face = m_physicalBCMap[bcId]->face;
3413 const MInt normalDir = face / 2;
3414 const MInt firstTangentialDir = (normalDir + 1) % nDim;
3415 const MInt normalDirStart = start[normalDir];
3416 const MInt firstTangentialStart = start[firstTangentialDir];
3417 const MInt firstTangentialEnd = end[firstTangentialDir];
3418 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
3419 // determine indices for direction help
3420 const MInt n = (face % 2) * 2 - 1; //-1,+1
3421 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
3422 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
3423 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
3424 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
3425
3426 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
3427 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
3428 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
3429 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
3430 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
3431
3432 const MFloat dxidx = m_cells->surfaceMetrics[normalDir * nDim + 0][cellIdA1];
3433 const MFloat dxidy = m_cells->surfaceMetrics[normalDir * nDim + 1][cellIdA1];
3434 // leaving domain of integration in positive coordinate direction,
3435 // therefore multiply with positive F1
3436 const MFloat gradxi = n * F1 / sqrt(dxidx * dxidx + dxidy * dxidy);
3437
3438 const MFloat rhoInner = m_cells->pvariables[PV->RHO][cellIdA1];
3439 const MFloat uInner = m_cells->pvariables[PV->U][cellIdA1];
3440 const MFloat vInner = m_cells->pvariables[PV->V][cellIdA1];
3441 const MFloat pInner = m_cells->pvariables[PV->P][cellIdA1];
3442
3443 const MFloat maContravariant = (dxidx * uInner + dxidy * vInner) * gradxi;
3444
3445
3446 if(maContravariant < F0) {
3447 // inflow
3448 const MFloat p = pInner;
3449 const MFloat rho = CV->rhoInfinity;
3450
3451 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3452 m_cells->pvariables[PV->U][cellIdG1] = PV->UInfinity;
3453 m_cells->pvariables[PV->V][cellIdG1] = PV->VInfinity;
3454 m_cells->pvariables[PV->P][cellIdG1] = p;
3455
3456 if(isRans) {
3457 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG1] = CV->ransInfinity[0] / CV->rhoInfinity;
3458 }
3459 } else {
3460 // outflow
3461 const MFloat p = PV->PInfinity;
3462 const MFloat rho = rhoInner;
3463
3464 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3465 m_cells->pvariables[PV->U][cellIdG1] = uInner;
3466 m_cells->pvariables[PV->V][cellIdG1] = vInner;
3467 m_cells->pvariables[PV->P][cellIdG1] = p;
3468 if(isRans) {
3469 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG1] =
3470 (F2 * m_cells->pvariables[PV->RANS_VAR[0]][cellIdA1] - m_cells->pvariables[PV->RANS_VAR[0]][cellIdA2]);
3471 }
3472 }
3473
3474 // extrapolate into second ghost cell
3475 for(MInt var = 0; var < PV->noVariables; var++) {
3476 m_cells->pvariables[var][cellIdG2] = F2 * m_cells->pvariables[var][cellIdG1] - m_cells->pvariables[var][cellIdA1];
3477 }
3478 }
3479}
StructuredCell::surfaceMetrics
MFloat ** surfaceMetrics
Definition: fvstructuredcell.h:44
maia::fc::cell::p
constexpr std::underlying_type< FcCell >::type p(const FcCell property)
Converts property name to underlying integer value.
Definition: fccellproperties.h:32

◆ bc2004()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2004 ( MInt  bcId)
virtual

Simplified characteristic approach (Whitfield) Still reflects pressure waves at the outflow rho=rho, u=u, v=v, w=w, p=p_inf

Author
Pascal Meysonnat

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3271 of file fvstructuredbndrycnd2d.cpp.

3271 {
3272 TRACE();
3273
3274 const MInt IJ[nDim] = {1, m_nCells[1]};
3275 const MFloat gamma = m_solver->m_gamma;
3276 MInt* start = m_physicalBCMap[bcId]->start1;
3277 MInt* end = m_physicalBCMap[bcId]->end1;
3278
3279 // Here we find out the normal direction of the
3280 // boundary and the tangential direction.
3281 // This way we can make a general formulation of
3282 // the boundary condition
3283 const MInt face = m_physicalBCMap[bcId]->face;
3284 const MInt normalDir = face / 2;
3285 const MInt firstTangentialDir = (normalDir + 1) % nDim;
3286 const MInt normalDirStart = start[normalDir];
3287 const MInt firstTangentialStart = start[firstTangentialDir];
3288 const MInt firstTangentialEnd = end[firstTangentialDir];
3289 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
3290 // determine indices for direction help
3291 const MInt n = (face % 2) * 2 - 1; //-1,+1
3292 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
3293 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
3294 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
3295 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
3296
3297 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
3298 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
3299 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
3300 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
3301 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
3302 const MFloat dxidx = m_cells->surfaceMetrics[normalDir * nDim + 0][cellIdA1];
3303 const MFloat dxidy = m_cells->surfaceMetrics[normalDir * nDim + 1][cellIdA1];
3304 // multiply with n, so it will be -1 or +1 depending if we enter
3305 // or leave the domain of integration in positive direction
3306 const MFloat gradxi = n * F1 / sqrt(dxidx * dxidx + dxidy * dxidy);
3307 const MFloat dxHelp = dxidx * gradxi;
3308 const MFloat dyHelp = dxidy * gradxi;
3309 // speed of sound
3310 const MFloat cBC = sqrt(gamma * m_cells->pvariables[PV->P][cellIdG1] / m_cells->pvariables[PV->RHO][cellIdG1]);
3311 const MFloat rhoBC = m_cells->pvariables[PV->RHO][cellIdG1];
3312 const MFloat rhoInner = m_cells->pvariables[PV->RHO][cellIdA1];
3313 const MFloat uInner = m_cells->pvariables[PV->U][cellIdA1];
3314 const MFloat vInner = m_cells->pvariables[PV->V][cellIdA1];
3315 const MFloat pInner = m_cells->pvariables[PV->P][cellIdA1];
3316 const MFloat maContravariant = (dxidx * uInner + dxidy * vInner - m_cells->dxt[normalDir][cellIdA1]) * gradxi;
3317 if(maContravariant < F0) {
3318 // inflow
3319 const MFloat p = F1B2
3320 * (pInner + PV->PInfinity
3321 + rhoBC * cBC * (dxHelp * (uInner - PV->UInfinity) + dyHelp * (vInner - PV->VInfinity)));
3322
3323 const MFloat rho = CV->rhoInfinity + (p - PV->PInfinity) / POW2(cBC);
3324 const MFloat help = (p - PV->PInfinity) / (rhoBC * cBC);
3325
3326 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3327 m_cells->pvariables[PV->U][cellIdG1] = (PV->UInfinity + help * dxHelp);
3328 m_cells->pvariables[PV->V][cellIdG1] = (PV->VInfinity + help * dyHelp);
3329 m_cells->pvariables[PV->P][cellIdG1] = p;
3330 if(isRans) {
3331 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3332 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdG1] = PV->ransInfinity[ransVarId];
3333 }
3334 }
3335 } else {
3336 // outflow
3337 const MFloat p = PV->PInfinity;
3338 const MFloat rho = rhoInner + (p - pInner) / POW2(cBC);
3339 const MFloat help = (p - pInner) / (rhoBC * cBC);
3340
3341 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3342 m_cells->pvariables[PV->U][cellIdG1] = (uInner - help * dxHelp);
3343 m_cells->pvariables[PV->V][cellIdG1] = (vInner - help * dyHelp);
3344 m_cells->pvariables[PV->P][cellIdG1] = p;
3345 if(isRans) {
3346 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3347 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdG1] =
3348 (F2 * m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdA1]
3349 - m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdA2]);
3350 }
3351 }
3352 }
3353
3354 // extrapolate into second ghost cell
3355 for(MInt var = 0; var < PV->noVariables; var++) {
3356 m_cells->pvariables[var][cellIdG2] = F2 * m_cells->pvariables[var][cellIdG1] - m_cells->pvariables[var][cellIdA1];
3357 }
3358 }
3359}
POW2
constexpr Real POW2(const Real x)
Definition: functions.h:119

◆ bc2005()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2005 ( MInt  bcId)
virtual

rho=rho, u=u, v=v, w=w, p=p

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3487 of file fvstructuredbndrycnd2d.cpp.

3487 {
3488 MInt* start = m_physicalBCMap[bcId]->start1;
3489 MInt* end = m_physicalBCMap[bcId]->end1;
3490 switch(m_physicalBCMap[bcId]->face) {
3491 case 1: {
3492 for(MInt j = start[1]; j < end[1]; j++) {
3493 for(MInt i = start[0]; i < end[0]; i++) {
3494 MInt cellId = cellIndex(i, j);
3495 MInt cellIdadj = cellIndex(i - 1, j);
3496
3497 for(MInt var = 0; var < PV->noVariables; var++) {
3498 m_cells->pvariables[var][cellId] = m_cells->pvariables[var][cellIdadj];
3499 }
3500 }
3501 }
3502 break;
3503 }
3504 case 3: {
3505 for(MInt j = start[1]; j < end[1]; j++) {
3506 for(MInt i = start[0]; i < end[0]; i++) {
3507 MInt cellId = cellIndex(i, j);
3508 MInt cellIdadj = cellIndex(i, j - 1);
3509 for(MInt var = 0; var < PV->noVariables; var++) {
3510 m_cells->pvariables[var][cellId] = m_cells->pvariables[var][cellIdadj];
3511 }
3512 }
3513 }
3514 break;
3515 }
3516 default: {
3517 mTerm(1, AT_, "Face direction not implemented");
3518 }
3519 }
3520}

◆ bc2006()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2006 ( MInt  bcId)
virtual

rho=rho, u=0, v=0, p=p_inf

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3528 of file fvstructuredbndrycnd2d.cpp.

3528 {
3529 const MInt IJ[nDim] = {1, m_nCells[1]};
3530 const MFloat gamma = m_solver->m_gamma;
3531 MInt* start = m_physicalBCMap[bcId]->start1;
3532 MInt* end = m_physicalBCMap[bcId]->end1;
3533
3534 // Here we find out the normal direction of the
3535 // boundary and the two tangential directions.
3536 // This way we can make a general formulation of
3537 // the boundary condition
3538 const MInt face = m_physicalBCMap[bcId]->face;
3539 const MInt normalDir = face / 2;
3540 const MInt firstTangentialDir = (normalDir + 1) % nDim;
3541 const MInt normalDirStart = start[normalDir];
3542 const MInt firstTangentialStart = start[firstTangentialDir];
3543 const MInt firstTangentialEnd = end[firstTangentialDir];
3544 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
3545
3546 const MInt n = (face % 2) * 2 - 1; //-1,+1
3547 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
3548 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
3549 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
3550 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
3551
3552 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
3553 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
3554 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
3555 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
3556 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
3557 const MFloat dxidx = m_cells->surfaceMetrics[normalDir * nDim + 0][cellIdA1];
3558 const MFloat dxidy = m_cells->surfaceMetrics[normalDir * nDim + 1][cellIdA1];
3559 // multiply with n, so it will be -1 or +1 depending if we enter
3560 // or leave the domain of integration in positive direction
3561 const MFloat gradxi = n * F1 / sqrt(dxidx * dxidx + dxidy * dxidy);
3562 const MFloat dxHelp = dxidx * gradxi;
3563 const MFloat dyHelp = dxidy * gradxi;
3564
3565 const MFloat cBC = sqrt(gamma * m_cells->pvariables[PV->P][cellIdG1] / m_cells->pvariables[PV->RHO][cellIdG1]);
3566 const MFloat rhoBC = m_cells->pvariables[PV->RHO][cellIdG1];
3567 const MFloat rhoInner = m_cells->pvariables[PV->RHO][cellIdA1];
3568 const MFloat uInner = m_cells->pvariables[PV->U][cellIdA1];
3569 const MFloat vInner = m_cells->pvariables[PV->V][cellIdA1];
3570 const MFloat pInner = m_cells->pvariables[PV->P][cellIdA1];
3571
3572 const MFloat maContravariant = (dxidx * uInner + dxidy * vInner - m_cells->dxt[normalDir][cellIdA1]) * gradxi;
3573
3574 if(maContravariant < F0) {
3575 // inflow
3576 const MFloat p =
3577 F1B2 * (pInner + PV->PInfinity + rhoBC * cBC * (dxHelp * (uInner - 0.0) + dyHelp * (vInner - 0.0)));
3578
3579 const MFloat rho = CV->rhoInfinity + (p - PV->PInfinity) / POW2(cBC);
3580 const MFloat help = (p - PV->PInfinity) / (rhoBC * cBC);
3581
3582 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3583 m_cells->pvariables[PV->U][cellIdG1] = (0.0 + help * dxHelp);
3584 m_cells->pvariables[PV->V][cellIdG1] = (0.0 + help * dyHelp);
3585 m_cells->pvariables[PV->P][cellIdG1] = p;
3586 if(isRans) {
3587 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG1] = PV->ransInfinity[0];
3588 }
3589 } else {
3590 // outflow
3591 const MFloat p = PV->PInfinity;
3592 const MFloat rho = rhoInner + (p - pInner) / POW2(cBC);
3593 const MFloat help = (p - pInner) / (rhoBC * cBC);
3594
3595 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3596 m_cells->pvariables[PV->U][cellIdG1] = (uInner - help * dxHelp);
3597 m_cells->pvariables[PV->V][cellIdG1] = (vInner - help * dyHelp);
3598 m_cells->pvariables[PV->P][cellIdG1] = p;
3599 if(isRans) {
3600 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG1] =
3601 (F2 * m_cells->pvariables[PV->RANS_VAR[0]][cellIdA1] - m_cells->pvariables[PV->RANS_VAR[0]][cellIdA2]);
3602 }
3603 }
3604
3605 // extrapolate into second ghost cell
3606 for(MInt var = 0; var < PV->noVariables; var++) {
3607 m_cells->pvariables[var][cellIdG2] = F2 * m_cells->pvariables[var][cellIdG1] - m_cells->pvariables[var][cellIdA1];
3608 }
3609 }
3610}

◆ bc2007()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2007 ( MInt  bcId)
virtual

rho=rho, u=u, v=v, w=w, p=p_inf

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3619 of file fvstructuredbndrycnd2d.cpp.

3619 {
3620 MInt* start = m_physicalBCMap[bcId]->start1;
3621 MInt* end = m_physicalBCMap[bcId]->end1;
3622
3623 switch(m_physicalBCMap[bcId]->face) {
3624 case 1: {
3625 for(MInt j = start[1]; j < end[1]; j++) {
3626 MInt cellId = cellIndex(start[0], j);
3627 MInt cellIdadj = cellIndex(start[0] - 1, j);
3628
3629 m_cells->pvariables[PV->RHO][cellId] = m_cells->pvariables[PV->RHO][cellIdadj];
3630 m_cells->pvariables[PV->U][cellId] = m_cells->pvariables[PV->U][cellIdadj];
3631 m_cells->pvariables[PV->V][cellId] = m_cells->pvariables[PV->V][cellIdadj];
3632 m_cells->pvariables[PV->P][cellId] = PV->PInfinity;
3633
3634 if(isRans) {
3635 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3636 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] =
3637 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdadj];
3638 }
3639 }
3640
3641 for(MInt var = 0; var < PV->noVariables; var++) {
3642 MInt cellIdP1 = cellIndex(start[0] + 1, j);
3643 m_cells->pvariables[var][cellIdP1] =
3644 2.0 * m_cells->pvariables[var][cellId] - m_cells->pvariables[var][cellIdadj];
3645 }
3646 }
3647 break;
3648 }
3649 case 3: {
3650 for(MInt i = start[0]; i < end[0]; i++) {
3651 MInt cellId = cellIndex(i, start[1]);
3652 MInt cellIdadj = cellIndex(i, start[1] - 1);
3653
3654 m_cells->pvariables[PV->RHO][cellId] = m_cells->pvariables[PV->RHO][cellIdadj];
3655 m_cells->pvariables[PV->U][cellId] = m_cells->pvariables[PV->U][cellIdadj];
3656 m_cells->pvariables[PV->V][cellId] = m_cells->pvariables[PV->V][cellIdadj];
3657 m_cells->pvariables[PV->P][cellId] = PV->PInfinity;
3658
3659 if(isRans) {
3660 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3661 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] =
3662 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdadj];
3663 }
3664 }
3665
3666 for(MInt var = 0; var < PV->noVariables; var++) {
3667 MInt cellIdP1 = cellIndex(i, start[1] + 1);
3668 m_cells->pvariables[var][cellIdP1] =
3669 2.0 * m_cells->pvariables[var][cellId] - m_cells->pvariables[var][cellIdadj];
3670 }
3671 }
3672
3673 break;
3674 }
3675
3676 default: {
3677 mTerm(1, AT_, "Face direction not implemented");
3678 }
3679 }
3680}

◆ bc2021()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2021 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3759 of file fvstructuredbndrycnd2d.cpp.

3759 {
3760 TRACE();
3761
3762 MInt* start = m_physicalBCMap[bcId]->start1;
3763 MInt* end = m_physicalBCMap[bcId]->end1;
3764 switch(m_physicalBCMap[bcId]->face) {
3765 case 0: {
3766 for(MInt j = start[1]; j < end[1]; j++) {
3767 const MInt cellIdG1 = cellIndex(1, j);
3768 const MInt cellIdG2 = cellIndex(0, j);
3769 const MInt cellIdA1 = cellIndex(2, j);
3770 // const MInt cellIdA2 = cellIndex(3,j);
3771
3772 const MFloat dxidx = m_cells->surfaceMetrics[0][cellIdA1];
3773 const MFloat dxidy = m_cells->surfaceMetrics[1][cellIdA1];
3774
3775 // multiply with n, so it will be -1 or +1 depending if we enter
3776 // or leave the domain of integration in positive direction
3777 const MFloat gradxi = -1 * F1 / sqrt(dxidx * dxidx + dxidy * dxidy);
3778
3779 const MFloat dxHelp = dxidx * gradxi;
3780 const MFloat dyHelp = dxidy * gradxi;
3781
3782 const MFloat cBC = sqrt(m_solver->m_gamma * pressure(cellIdG1) / m_cells->pvariables[PV->RHO][cellIdG1]);
3783 const MFloat rhoBC = m_cells->pvariables[PV->RHO][cellIdG1];
3784
3785 const MFloat uInner = m_cells->pvariables[PV->U][cellIdA1];
3786 const MFloat vInner = m_cells->pvariables[PV->V][cellIdA1];
3787 const MFloat pInner = pressure(cellIdA1);
3788
3789 const MFloat uInflow = PV->UInfinity * (m_cells->coordinates[1][cellIdG1] * m_bc2021Gradient);
3790 const MFloat vInflow = 0.0;
3791
3792 // inflow
3793 const MFloat p =
3794 F1B2 * (pInner + PV->PInfinity + rhoBC * cBC * (dxHelp * (uInner - uInflow) + dyHelp * (vInner - vInflow)));
3795
3796 const MFloat rho = CV->rhoInfinity + (p - PV->PInfinity) / POW2(cBC);
3797 const MFloat help = (p - PV->PInfinity) / (rhoBC * cBC);
3798
3799 m_cells->pvariables[PV->RHO][cellIdG1] = rho;
3800 m_cells->pvariables[PV->U][cellIdG1] = uInflow + help * dxHelp;
3801 m_cells->pvariables[PV->V][cellIdG1] = vInflow + help * dyHelp;
3802 m_cells->pvariables[PV->P][cellIdG1] = p;
3803
3804 if(isRans) {
3805 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3806 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdG1] = PV->ransInfinity[ransVarId];
3807 }
3808 }
3809
3810 // extrapolate into second ghost cell
3811 for(MInt var = 0; var < PV->noVariables; var++) {
3812 m_cells->pvariables[var][cellIdG2] =
3813 F2 * m_cells->pvariables[var][cellIdG1] - m_cells->pvariables[var][cellIdA1];
3814 }
3815 }
3816 break;
3817 }
3818 default:
3819
3820 {
3821 mTerm(1, AT_, "Face direction not implemented)");
3822 }
3823 }
3824}
StructuredCell::coordinates
MFloat ** coordinates
Definition: fvstructuredcell.h:15

◆ bc2199()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2199 ( MInt  bcId)
virtual

!!!!!!!!! check for the gamma stuff in tfs

!!!!!!!check the sign convention from tfs

!!!!!! check for ii1

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3827 of file fvstructuredbndrycnd2d.cpp.

3827 {
3828 TRACE();
3830 const MFloat gamma = m_solver->m_gamma;
3831 const MFloat gammaMinusOne = m_solver->m_gamma - F1;
3832 const MFloat fgamma = F1 / gamma;
3833 MInt* start = m_physicalBCMap[bcId]->start1;
3834 MInt* end = m_physicalBCMap[bcId]->end1;
3835 switch(m_physicalBCMap[bcId]->face) {
3836 case 0: {
3837 MFloat mflux = F0;
3838 MFloat surface = F0;
3839 for(MInt j = start[1]; j < end[1]; j++) {
3840 const MInt cellId = cellIndex(m_noGhostLayers, j); // first inner line
3841 const MInt pointId = getPointIdFromCell(m_noGhostLayers, j);
3842 const MInt pointIdP1 = getPointIdFromPoint(pointId, 0, 1);
3844 const MFloat dx = m_grid->m_coordinates[0][pointIdP1] - m_grid->m_coordinates[0][pointId];
3845 const MFloat dy = m_grid->m_coordinates[1][pointIdP1] - m_grid->m_coordinates[1][pointId];
3846 const MFloat rhou = m_cells->pvariables[PV->RHO][cellId] * m_cells->pvariables[PV->U][cellId];
3847 const MFloat rhov = m_cells->pvariables[PV->RHO][cellId] * m_cells->pvariables[PV->V][cellId];
3848 mflux += dx * rhou + dy * rhov;
3849 surface += sqrt(POW2(dx) + POW2(dy));
3850 }
3851 // get cellId at the domain top
3853 MInt cellIdII1 = cellIndex(m_noGhostLayers, end[1] - m_noGhostLayers);
3854 // get averaged flux
3855 const MFloat rhouAVG = mflux / surface;
3856 // determine the pressure at the top at the iner line
3857 MFloat pressure = min(F1, max(F0, m_cells->pvariables[PV->P][cellIdII1] * gamma));
3858 // fix poin iteration apparently converges fast
3859 // from schlichting and trockenbrot, page 155
3860 // ru=sqrt(2.*fgamm1)*pp**fgam*sqrt(f1-pp**(gamm1*fgam))
3861 for(MInt i = 0; i < 20; i++) {
3862 pressure = pow((F1 - pow((rhouAVG / (sqrt(F2 * gammaMinusOne) * pow(pressure, fgamma))), F2)),
3863 (gammaMinusOne * fgamma));
3864 }
3865 pressure *= fgamma;
3866
3867 // scale the velocity profile
3868
3869 break;
3870 }
3871 default: {
3872 mTerm(1, AT_, "Face direction not implemented)");
3873 }
3874 }
3875}

◆ bc2402()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2402 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3879 of file fvstructuredbndrycnd2d.cpp.

3879 {
3880 TRACE();
3881
3882 MInt* start = m_physicalBCMap[bcId]->start1;
3883 MInt* end = m_physicalBCMap[bcId]->end1;
3884
3885 // first we need the average pressure and Temperature
3886 //==> calculate from the field
3887 //==> average over the surface
3888 //==> store in variable
3889
3890 MInt Id = m_channelSurfaceIndexMap[bcId];
3891 MInt* startface = &m_solver->m_windowInfo->channelSurfaceIndices[Id]->start1[0];
3892 MInt* endface = &m_solver->m_windowInfo->channelSurfaceIndices[Id]->end1[0];
3893 // global Pressure and Temperature at in-/outlet
3894 MFloat globalTin[2] = {F0, F0};
3895 MFloat globalTout[2] = {F0, F0};
3896 MFloat globalPin[2] = {F0, F0};
3897 MFloat globalPout[2] = {F0, F0};
3898 cout.precision(10);
3899
3900 // find out which face it is
3901 switch(m_physicalBCMap[bcId]->face) {
3902 case 0: {
3903 // average the pressure
3904 // use values from the inside to determine the pressure at the face!!!
3905 MFloat surface = F0;
3906 MFloat localPin[2] = {F0, F0};
3907 MFloat localTin[2] = {F0, F0};
3908 MFloat localVel = F0, globalVel = F0;
3909 MFloat localMassFlux = F0, globalMassFlux = F0;
3910
3911 for(MInt j = startface[1]; j < endface[1]; j++) {
3912 MInt ii = end[0] - 1;
3913 MInt cellIdP1 = cellIndex(ii + 1, j);
3914 MInt cellIdP2 = cellIndex(ii + 2, j);
3915
3916 surface = sqrt(POW2(m_cells->surfaceMetrics[0][cellIdP1]) + POW2(m_cells->surfaceMetrics[1][cellIdP1]));
3917 localPin[0] += surface * m_cells->pvariables[PV->P][cellIdP1];
3918 localPin[1] += surface * m_cells->pvariables[PV->P][cellIdP2];
3919 localTin[0] += surface * temperature(cellIdP1);
3920 localTin[1] += surface * temperature(cellIdP2);
3921 localVel += surface * (m_cells->pvariables[PV->U][cellIdP1]);
3922 localMassFlux += surface * (m_cells->pvariables[PV->U][cellIdP1] * m_cells->pvariables[PV->RHO][cellIdP1]);
3923 }
3924
3925 // next now that the pressure and temperature are known:
3926 // make it a global variable for the complete inflow plane
3927 MPI_Allreduce(localPin, globalPin, 2, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelIn, AT_, "localPin",
3928 "globalPin");
3929 MPI_Allreduce(localTin, globalTin, 2, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelIn, AT_, "localTin",
3930 "globalTin");
3931 MPI_Allreduce(&localMassFlux, &globalMassFlux, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelIn, AT_,
3932 "localMassFlux", "globalMassFlux");
3933 MPI_Allreduce(&localVel, &globalVel, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelIn, AT_, "localVel",
3934 "globalVel");
3935
3936 MPI_Barrier(m_solver->m_commChannelIn, AT_);
3937 globalPin[0] /= m_channelSurfaceIn;
3938 globalTin[0] /= m_channelSurfaceIn;
3939 globalPin[1] /= m_channelSurfaceIn;
3940 globalTin[1] /= m_channelSurfaceIn;
3941 globalMassFlux /= m_channelSurfaceIn;
3942 globalVel /= m_channelSurfaceIn;
3943 MFloat currentRe = m_solver->m_Re / PV->UInfinity * globalVel;
3944
3945 // set a Barrier to ensure that all domains having a channel side are at same time level
3946 // now make the variables known everywhere in the channel world
3947 MPI_Barrier(m_solver->m_commChannelWorld, AT_);
3948 MPI_Bcast(globalTin, 2, MPI_DOUBLE, m_solver->m_channelRoots[2], m_solver->m_commChannelWorld, AT_, "globalTin");
3949 MPI_Bcast(globalPin, 2, MPI_DOUBLE, m_solver->m_channelRoots[2], m_solver->m_commChannelWorld, AT_, "globalPin");
3950 MPI_Bcast(globalTout, 2, MPI_DOUBLE, m_solver->m_channelRoots[3], m_solver->m_commChannelWorld, AT_,
3951 "globalTout");
3952 MPI_Bcast(globalPout, 2, MPI_DOUBLE, m_solver->m_channelRoots[3], m_solver->m_commChannelWorld, AT_,
3953 "globalPout");
3954 // now every value is known and can be used to apply the BC !!!!
3955 // cannot take values from above at here the window is extended to the ghost layers
3956
3957 if(globalTimeStep > 1 && globalTimeStep % 50 == 0) {
3958 if(m_channelInflowRank == 0 && globalTimeStep % m_solver->m_residualInterval == 0 && m_solver->m_RKStep == 0) {
3959 cout.precision(6);
3960
3961 FILE* f_channel;
3962 f_channel = fopen("./massflow.dat", "a+");
3963 fprintf(f_channel, "%d", globalTimeStep);
3964 fprintf(f_channel, " %f", m_solver->m_physicalTime);
3965 fprintf(f_channel, " %f", m_solver->m_time);
3966 fprintf(f_channel, " %f", m_solver->m_timeStep);
3967 fprintf(f_channel, " %f", currentRe);
3968 fprintf(f_channel, " %f", globalMassFlux);
3969 fprintf(f_channel, " %f", globalPin[0]);
3970 fprintf(f_channel, " %f", globalPin[1]);
3971 fprintf(f_channel, " %f", globalTin[0]);
3972 fprintf(f_channel, " %f", globalTout[0]);
3973 fprintf(f_channel, " %f", globalPout[0]);
3974 fprintf(f_channel, " %f", globalPout[1]);
3975 fprintf(f_channel, " %f", (globalTin[0] - globalTout[1]));
3976 fprintf(f_channel, "\n");
3977 fclose(f_channel);
3978 }
3979 }
3980
3981 // If fully periodic, do nothing
3982 if(m_solver->m_channelFullyPeriodic) {
3983 m_solver->m_inflowVelAvg = globalVel;
3984 return;
3985 }
3986
3987 for(MInt j = start[1]; j < end[1]; j++) {
3988 MInt i = 1;
3989 const MInt cellId = cellIndex(i, j);
3990
3991 MFloat pressureFluctuation = m_cells->pvariables[PV->P][cellId] - globalPout[1];
3992 MFloat x = m_cells->coordinates[0][cellId];
3993 MFloat pressureInflowMean =
3994 m_solver->m_deltaP / m_solver->m_channelLength * (x - m_solver->m_channelInflowPlaneCoordinate)
3995 + PV->PInfinity;
3996 MFloat pressureNew = pressureInflowMean + pressureFluctuation;
3997 MFloat temperatureFlucOutflow = temperature(cellId) - globalTout[1];
3998 MFloat temperatureNew = PV->TInfinity + temperatureFlucOutflow;
3999 MFloat rhoNew = m_solver->m_gamma * pressureNew / temperatureNew;
4000
4001 m_cells->pvariables[PV->RHO][cellId] = rhoNew;
4002 m_cells->pvariables[PV->P][cellId] = pressureNew;
4003
4004 for(MInt var = 0; var < PV->noVariables; var++) {
4005 m_cells->pvariables[var][cellIndex(start[0], j)] = 2.0 * m_cells->pvariables[var][cellIndex(start[0] + 1, j)]
4006 - m_cells->pvariables[var][cellIndex(start[0] + 2, j)];
4007 }
4008 }
4009 break;
4010 }
4011 case 1: {
4012 // average the pressure
4013 // use values from the inside to determine the pressure at the face!!!
4014 MFloat surface = F0;
4015 MFloat localPout[] = {F0, F0};
4016 MFloat localTout[] = {F0, F0};
4017 for(MInt j = startface[1]; j < endface[1]; j++) {
4018 MInt ii = start[0];
4019 MInt cellIdM1 = cellIndex(ii - 1, j);
4020 MInt cellIdM2 = cellIndex(ii - 2, j);
4021 surface = sqrt(POW2(m_cells->surfaceMetrics[0][cellIdM1]) + POW2(m_cells->surfaceMetrics[1][cellIdM1]));
4022 localPout[0] += surface * m_cells->pvariables[PV->P][cellIdM2];
4023 localPout[1] += surface * m_cells->pvariables[PV->P][cellIdM1];
4024 localTout[0] += surface * temperature(cellIdM2);
4025 localTout[1] += surface * temperature(cellIdM1);
4026 }
4027
4028 // next now that the pressure and temperature are known:
4029 // make it a global variable for the complete outflow plane
4030 MPI_Allreduce(localPout, globalPout, 2, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelOut, AT_, "localPout",
4031 "globalPout");
4032 MPI_Allreduce(localTout, globalTout, 2, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelOut, AT_, "localTout",
4033 "globalTout");
4034
4035 MPI_Barrier(m_solver->m_commChannelOut, AT_);
4036 globalPout[0] /= m_channelSurfaceOut;
4037 globalTout[0] /= m_channelSurfaceOut;
4038 globalPout[1] /= m_channelSurfaceOut;
4039 globalTout[1] /= m_channelSurfaceOut;
4040 // set a Barrier to ensure that all domains having a channel side are at same time level
4041 // now make the variables known everywhere in the channel world
4042 MPI_Barrier(m_solver->m_commChannelWorld, AT_);
4043 MPI_Bcast(globalTin, 2, MPI_DOUBLE, m_solver->m_channelRoots[2], m_solver->m_commChannelWorld, AT_, "globalTin");
4044 MPI_Bcast(globalPin, 2, MPI_DOUBLE, m_solver->m_channelRoots[2], m_solver->m_commChannelWorld, AT_, "globalPin");
4045 MPI_Bcast(globalTout, 2, MPI_DOUBLE, m_solver->m_channelRoots[3], m_solver->m_commChannelWorld, AT_,
4046 "globalTout");
4047 MPI_Bcast(globalPout, 2, MPI_DOUBLE, m_solver->m_channelRoots[3], m_solver->m_commChannelWorld, AT_,
4048 "globalPout");
4049 // now every value is known and can be used to apply the BC !!!!
4050 // cannot take values from above at here the window is extended to the ghost layers
4051
4052 // If fully periodic, do nothing; return after mpi-stuff to avoid blocking
4053 if(m_solver->m_channelFullyPeriodic) {
4054 return;
4055 }
4056
4057 for(MInt j = start[1]; j < end[1]; j++) {
4058 MInt i = start[0];
4059 const MInt cellId = cellIndex(i, j);
4060
4061 MFloat pressureFluctuation = m_cells->pvariables[PV->P][cellId] - globalPin[0];
4062 MFloat x = m_cells->coordinates[0][cellId];
4063 MFloat pressureOutflowMean =
4064 m_solver->m_deltaP / m_solver->m_channelLength * (x - m_solver->m_channelInflowPlaneCoordinate)
4065 + PV->PInfinity;
4066 MFloat pressureNew = pressureOutflowMean + pressureFluctuation;
4067 MFloat deltaT = (globalTin[0] - globalTout[1]);
4068 MFloat temperatureInflow = temperature(cellId);
4069 MFloat temperatureNew = temperatureInflow - deltaT;
4070 MFloat rhoNew = m_solver->m_gamma * pressureNew / temperatureNew;
4071
4072 m_cells->pvariables[PV->RHO][cellId] = rhoNew;
4073 m_cells->pvariables[PV->P][cellId] = pressureNew;
4074
4075 for(MInt var = 0; var < PV->noVariables; var++) {
4076 m_cells->pvariables[var][cellIndex(start[0] + 1, j)] = 2.0 * m_cells->pvariables[var][cellIndex(start[0], j)]
4077 - m_cells->pvariables[var][cellIndex(start[0] - 1, j)];
4078 }
4079 }
4080 break;
4081 }
4082 default: {
4083 mTerm(1, AT_, "Face not implemented");
4084 break;
4085 }
4086 }
4087}
FvStructuredSolver::m_windowInfo
std::unique_ptr< FvStructuredSolverWindowInfo< nDim > > m_windowInfo
Definition: fvstructuredsolver.h:419
Solver::m_residualInterval
MInt m_residualInterval
The number of timesteps before writing the next residual.
Definition: solver.h:87
Solver::m_Re
MFloat m_Re
the Reynolds number
Definition: solver.h:68
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
MPI_Barrier
int MPI_Barrier(MPI_Comm comm, const MString &name)
same as MPI_Barrier
Definition: mpioverride.cpp:874
MPI_Allreduce
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
Definition: mpioverride.cpp:952
MPI_Bcast
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
Definition: mpioverride.cpp:1114

◆ bc2510()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2510 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 4092 of file fvstructuredbndrycnd2d.cpp.

4092 {
4093 (void)bcId;
4094
4095 cout.precision(8);
4096 // MInt* start = m_physicalBCMap[bcId]->start1;
4097 // MInt* end = m_physicalBCMap[bcId]->end1;
4098 // Communication between the Recycling and Inlet station.
4099 //___COMMGr(oup)______________________
4100 //|______________|_________________| |
4101 //|Inlet station |Recycling station| |
4102 //| | | |
4103 //|______________|_________________| |
4104 //|__________________________________|
4105 //
4106 // infographic of the Groups
4107
4108 // MFloat delta_inMax = delta_in+2.5*delta_in; //limit the integration height
4109 const MFloat gamma = m_solver->m_gamma;
4110 const MFloat gammaMinusOne = gamma - F1;
4111
4112 const MFloat rescalEPS = pow(10, -16.0);
4113 const MFloat alpha = 4.0;
4114 const MFloat b = 0.2;
4115 const MFloat rc = pow(m_solver->m_Pr, F1B3);
4116 const MFloat ctema = F1B2 * gammaMinusOne * POW2(m_solver->m_Ma) * rc;
4117 const MFloat maxIntegrationHeight = 2.0 * m_rescalingBLT;
4118
4119 // van Driest constant & transformed velocity
4120 const MFloat b_vd = sqrt(ctema / (F1 + ctema));
4121 const MFloat uvd8 = PV->UInfinity * asin(b_vd) / b_vd;
4122 const MInt i = m_noGhostLayers - 1;
4123
4127
4128 // allocate space in k direction (all k-Cells)
4129 MFloatScratchSpace thetaLocal(2, AT_, "thetaLocalIn");
4130 thetaLocal.fill(F0); // initialize scratch space
4131
4132 MFloatScratchSpace thetaGlobal(2, AT_, "thetaGlobalIn");
4133 thetaGlobal.fill(F0);
4134
4135 // compute the local moment thickness j=direction of integration
4136 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; ++j) {
4137 const MInt cellId = cellIndex(i, j);
4138 const MInt pointIdM1 = getPointIdFromCell(i, j);
4139 const MInt pointIdP1 = getPointIdFromPoint(pointIdM1, 0, 1);
4140
4141 if(m_grid->m_coordinates[1][pointIdM1] > maxIntegrationHeight) {
4142 continue;
4143 }
4144
4145 const MFloat urat = m_cells->pvariables[PV->U][cellId] / PV->UInfinity;
4146 const MFloat momThick = m_cells->pvariables[PV->U][cellId] * m_cells->pvariables[PV->RHO][cellId] * fabs(F1 - urat)
4147 / (CV->rhoUInfinity);
4148 // integrate normal to the wall
4149 const MFloat ydist = m_grid->m_coordinates[1][pointIdP1] - m_grid->m_coordinates[1][pointIdM1];
4150 thetaLocal(0) += momThick * ydist;
4151 }
4152
4153 MPI_Allreduce(thetaLocal.begin(), thetaGlobal.begin(), 2, MPI_DOUBLE, MPI_SUM, m_solver->m_rescalingCommGrComm, AT_,
4154 "thetaLocal.begin()", "thetaGlobal.begin()");
4155
4156 thetaGlobal(0) = mMax(thetaGlobal(0), 0.0000001);
4157 thetaGlobal(1) = mMax(thetaGlobal(1), 0.0000001);
4158
4159 if(globalTimeStep % 50 == 0 && m_solver->m_RKStep == 0
4160 && m_solver->domainId() == m_solver->m_rescalingCommGrRootGlobal) {
4161 cout << m_solver->domainId() << " ThetaInflow " << thetaGlobal(0) << " ThetaRecyclingStation " << thetaGlobal(1)
4162 << endl;
4163
4164 FILE* f_channel;
4165 f_channel = fopen("./theta_inflow.dat", "a+");
4166 fprintf(f_channel, "%d", globalTimeStep);
4167 fprintf(f_channel, " %f", m_solver->m_physicalTime);
4168 fprintf(f_channel, " %f", m_solver->m_time);
4169 fprintf(f_channel, " %f", m_solver->m_timeStep);
4170 fprintf(f_channel, " %f", thetaGlobal[0]);
4171 fprintf(f_channel, " %f", thetaGlobal[1]);
4172 fprintf(f_channel, "\n");
4173 fclose(f_channel);
4174 }
4175
4179 const MInt noWallProperties = 3;
4180 MFloatScratchSpace wallPropertiesLocal(noWallProperties, AT_, "wallPropertiesLocalInlet");
4181 MFloatScratchSpace wallProperties(noWallProperties, AT_, "wallPropertiesInlet");
4182 wallPropertiesLocal.fill(F0);
4183 wallProperties.fill(F0);
4184
4185 MPI_Allreduce(wallPropertiesLocal.begin(), wallProperties.begin(), noWallProperties, MPI_DOUBLE, MPI_SUM,
4186 m_solver->m_rescalingCommGrComm, AT_, "wallPropertiesLocal.begin()", "wallProperties.begin()");
4187
4191
4192 MFloat utauIn = F0;
4193 MFloat gams = F0;
4194
4195 MFloat utauRe = wallProperties(0);
4196 // estimate the friction velocity at the inlet
4197 // according to standard power law approximations
4198 // utau_in = utau_re*(theta_re/theta_in)**(1/2*(n-1))
4199 // where theta is the momentum thichness
4200 // see Thomas S. Lund, p241
4201
4202 // when take into account the variance of wall density
4203 // utau_in = utau_re*(rho_wall_re/rho_wall_in)**0.5
4204 //* (theta_re/theta_in)**(1/2*(n-1))
4205 // here n = 5
4206 const MFloat n = 5.0;
4207 const MFloat facc = F1 / (F2 * (n - F1));
4208
4209 gams = pow(thetaGlobal(1) / fabs(thetaGlobal(0)), facc);
4210 gams = mMin(gams, 1.8);
4211 utauIn = utauRe * gams;
4212
4213 MFloatScratchSpace coordInInner(m_nCells[0], AT_, "coordInInner");
4214 MFloatScratchSpace coordInOuter(m_nCells[0], AT_, "coordInOuter");
4215
4216 for(MInt j = 0; j < m_nCells[0]; ++j) {
4217 const MInt cellId = cellIndex(i, j);
4218 const MFloat rho = m_cells->pvariables[PV->RHO][cellId];
4219 const MFloat frho = F1 / rho;
4220 const MFloat p = m_cells->pvariables[PV->P][cellId];
4221 const MFloat temp = p * gamma * frho;
4222 const MFloat mu = SUTHERLANDLAW(temp);
4223
4224 coordInInner(j) = utauIn * rho * m_cells->coordinates[1][cellId] / (mu * sqrt(m_solver->m_Re0));
4225 coordInOuter(j) = m_cells->coordinates[1][cellId] * rho / (m_rescalingBLT * CV->rhoInfinity);
4226 }
4227
4228 const MInt wallLocalOffset = m_solver->m_nOffsetCells[0]; // Offset in j-direction
4229 MFloat tempWallInletLocal = F0;
4230 MFloat tempWallInletGlobal = F0;
4231 // determine the wall stuff if wall is contained whithin the partition
4232 if(wallLocalOffset == 0 && m_solver->m_nActiveCells[0] >= m_noGhostLayers) {
4233 const MInt cellId = cellIndex(i, m_noGhostLayers);
4234 tempWallInletLocal = temperature(cellId);
4235 }
4236
4237 MPI_Allreduce(&tempWallInletLocal, &tempWallInletGlobal, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_rescalingCommGrComm, AT_,
4238 "tempWallInletLocal", "tempWallInletGlobal");
4239
4243 const MInt noVariables = PV->noVariables + 1;
4244 MInt totalCells = m_grid->getMyBlockNoCells(0) + 1;
4245 MFloatScratchSpace varSliceLocal(noVariables, totalCells, AT_, "varSliceLocal");
4246 MFloatScratchSpace varSlice(noVariables, totalCells, AT_, "varSlice");
4247
4248 // we are at the inlet, only fill with zeros
4249 varSlice.fill(F0);
4250 varSliceLocal.fill(F0);
4251
4252 MPI_Allreduce(varSliceLocal.begin(), varSlice.begin(), noVariables * totalCells, MPI_DOUBLE, MPI_SUM,
4253 m_solver->m_rescalingCommGrComm, AT_, "varSliceLocal.begin()", "varSlice.begin()");
4254
4255 // if(globalTimeStep%5==0&&m_solver->m_RKStep==0&&m_solver->domainId() == m_solver->m_rescalingCommGrRootGlobal) {
4256 // cout << m_solver->domainId()<< " ThetaInflow " << thetaGlobal(0) <<" ThetaRecyclingStation " << thetaGlobal(1) <<
4257 // endl;
4258
4259 // FILE* f_channel;
4260 // f_channel = fopen("./varslice.dat", "w");
4261 // for(MInt jj=0; jj<totalCells-1; ++jj) {
4262 // for(MInt var=0; var<6; var++) {
4263 // fprintf(f_channel, " %f", varSlice(var, jj));
4264 // }
4265 // fprintf(f_channel, "\n");
4266 // }
4267 // fclose(f_channel);
4268 // }
4269
4273
4274 const MFloat ctem1 = (F1 + ctema) * (F1 - POW2(gams));
4275 const MFloat ctem2 = F2 * ctema * gams * (F1 - gams);
4276 const MFloat ctem3 = (F1 - gams) * (F1 + gams + F2 * ctema * gams);
4277
4278 MInt jStart = 0;
4279 MInt jEnd = m_nCells[0];
4280
4281 if(m_solver->m_nOffsetCells[0] == 0) {
4282 jStart = m_noGhostLayers;
4283 }
4284 if(m_solver->m_nOffsetCells[0] + m_solver->m_nActiveCells[0] == m_grid->getMyBlockNoCells(0)) {
4285 jEnd = m_nCells[0] - m_noGhostLayers;
4286 }
4287
4288 MFloat blEdgeVValueLocal = F0;
4289 MBool edgePointIsSet = false;
4290 MInt edgePointJ = 0;
4291
4292
4293 for(MInt j = jStart; j < jEnd; ++j) {
4294 const MInt cellId = cellIndex(i, j);
4295
4296 if(j > 0) {
4297 if(coordInOuter(j - 1) < 1.05 && coordInOuter(j) >= 1.05) {
4298 blEdgeVValueLocal = m_cells->pvariables[PV->V][cellIndex(i, j - 1)];
4299 }
4300 }
4301
4302 if(coordInOuter(j) < 1.05) {
4303 MFloat uInner = F0, vInner = F0, TInner = F0, mutInner = F0;
4304 MFloat uOuter = F0, vOuter = F0, TOuter = F0, mutOuter = F0;
4305 const MFloat count = alpha * (coordInOuter(j) - b);
4306 const MFloat denom = (F1 - F2 * b) * coordInOuter(j) + b;
4307 const MFloat ratio = count / denom;
4308 const MFloat wfun = F1B2 * (F1 + tanh(ratio) / tanh(alpha));
4309
4310 for(MInt jj = 0; jj < totalCells - 1; ++jj) {
4311 const MInt localId = jj;
4312 const MInt localIdP1 = jj + 1;
4313
4314 const MFloat yInnerRe = varSlice(3, localId);
4315 const MFloat yInnerReP1 = varSlice(3, localIdP1);
4316
4317 if((yInnerRe - coordInInner(j)) < rescalEPS && yInnerReP1 > coordInInner(j)) {
4318 const MFloat dy1 = coordInInner(j) - yInnerRe;
4319 const MFloat dy2 = yInnerReP1 - coordInInner(j);
4320 const MFloat dy = yInnerReP1 - yInnerRe;
4321
4322 const MFloat u = varSlice(0, localId);
4323 const MFloat uP1 = varSlice(0, localIdP1);
4324 const MFloat v = varSlice(1, localId);
4325 const MFloat vP1 = varSlice(1, localIdP1);
4326 const MFloat t = varSlice(2, localId);
4327 const MFloat tP1 = varSlice(2, localIdP1);
4328 const MFloat mut = varSlice(5, localId);
4329 const MFloat mutP1 = varSlice(5, localIdP1);
4330 uInner = (uP1 * dy1 + u * dy2) / dy;
4331 vInner = (vP1 * dy1 + v * dy2) / dy;
4332 TInner = (tP1 * dy1 + t * dy2) / dy;
4333 mutInner = (mutP1 * dy1 + mut * dy2) / dy;
4334 }
4335 }
4336
4337 // outer region
4338 for(MInt jj = 0; jj < totalCells - 1; ++jj) {
4339 const MInt localId = jj;
4340 const MInt localIdP1 = jj + 1;
4341
4342 const MFloat yOuterRe = varSlice(4, localId);
4343 const MFloat yOuterReP1 = varSlice(4, localIdP1);
4344
4345 if((yOuterRe - coordInOuter(j)) < rescalEPS && yOuterReP1 > coordInOuter(j)) {
4346 const MFloat dy1 = coordInOuter(j) - yOuterRe;
4347 const MFloat dy2 = yOuterReP1 - coordInOuter(j);
4348 const MFloat dy = yOuterReP1 - yOuterRe;
4349
4350 const MFloat u = varSlice(0, localId);
4351 const MFloat uP1 = varSlice(0, localIdP1);
4352 const MFloat v = varSlice(1, localId);
4353 const MFloat vP1 = varSlice(1, localIdP1);
4354 const MFloat t = varSlice(2, localId);
4355 const MFloat tP1 = varSlice(2, localIdP1);
4356 const MFloat mut = varSlice(5, localId);
4357 const MFloat mutP1 = varSlice(5, localIdP1);
4358 uOuter = (uP1 * dy1 + u * dy2) / dy;
4359 vOuter = (vP1 * dy1 + v * dy2) / dy;
4360 TOuter = (tP1 * dy1 + t * dy2) / dy;
4361 mutOuter = (mutP1 * dy1 + mut * dy2) / dy;
4362 }
4363 }
4364
4365 const MFloat TInnerA = POW2(gams) * TInner + ctem1 * PV->TInfinity;
4366 const MFloat TOuterA = POW2(gams) * TOuter - (ctem2 * (uOuter / PV->UInfinity) - ctem3) * PV->TInfinity;
4367
4368 // van Driest transformation
4369 const MFloat uvdInner = PV->UInfinity * asin(b_vd * uInner / PV->UInfinity) / b_vd;
4370 const MFloat uvdOuter = PV->UInfinity * asin(b_vd * uOuter / PV->UInfinity) / b_vd;
4371
4372 // scaling of transformed inner and outer velocities
4373 // use uvd8, the van Driest transformed u8 value
4374 uInner = gams * uvdInner;
4375 uOuter = gams * uvdOuter + (F1 - gams) * uvd8;
4376
4377 uInner = PV->UInfinity * sin(b_vd * uInner / PV->UInfinity) / b_vd;
4378 uOuter = PV->UInfinity * sin(b_vd * uOuter / PV->UInfinity) / b_vd;
4379
4380 // turbulent viscosity
4381 const MFloat tempWallInlet = tempWallInletGlobal;
4382 const MFloat tempWallRecycling = wallProperties(2);
4383
4384 const MFloat viscWallInlet = SUTHERLANDLAW(tempWallInlet);
4385 const MFloat viscWallRecycling = SUTHERLANDLAW(tempWallRecycling);
4386 const MFloat thetaInlet = thetaGlobal(0);
4387 const MFloat thetaRecycling = thetaGlobal(1);
4388 mutInner = mutInner * (viscWallInlet / viscWallRecycling);
4389 mutOuter = mutOuter * gams * (thetaInlet / thetaRecycling);
4390
4391 const MFloat uMean = uInner * (F1 - wfun) + uOuter * wfun;
4392 const MFloat vMean = vInner * (F1 - wfun) + vOuter * wfun;
4393 const MFloat tMean = TInnerA * (F1 - wfun) + TOuterA * wfun;
4394 MFloat mutMean = mutInner * (F1 - wfun) + mutOuter * wfun;
4395
4396 const MFloat clebf = 6.6;
4397 const MFloat blt = m_rescalingBLT;
4398 const MFloat cleb = F1 / (F1 + pow((m_cells->coordinates[1][cellId] / (clebf * blt)), 6.0));
4399
4400 const MFloat pres = PV->PInfinity;
4401 const MFloat rhoIn = gamma * pres / tMean;
4402
4403 m_cells->pvariables[PV->RHO][cellId] = rhoIn * cleb;
4404 m_cells->pvariables[PV->U][cellId] = uMean;
4405 m_cells->pvariables[PV->V][cellId] = vMean;
4406 m_cells->pvariables[PV->P][cellId] = pres;
4407 m_cells->pvariables[PV->RANS_VAR[0]][cellId] = mutMean / rhoIn;
4408
4409 } else {
4410 if(!edgePointIsSet) {
4411 edgePointJ = j;
4412 edgePointIsSet = true;
4413 }
4414 // const MFloat pres = pressure(cellIndex(m_noGhostLayers,j,k)); //for supersonic PV->PInfinity
4415 const MFloat pres = PV->PInfinity;
4416 const MFloat rhoIn = gamma * pres / PV->TInfinity;
4417
4418 const MFloat uMean = PV->UInfinity;
4419 const MFloat vMean = PV->VInfinity;
4420
4421 m_cells->pvariables[PV->RHO][cellId] = rhoIn;
4422 m_cells->pvariables[PV->U][cellId] = uMean;
4423 m_cells->pvariables[PV->V][cellId] = vMean; // m_cells->pvariables[PV->V][cellIndex(i,j-1)]*rhoIn;
4424 m_cells->pvariables[PV->P][cellId] = pres;
4425 m_cells->pvariables[PV->RANS_VAR[0]][cellId] = PV->ransInfinity[0];
4426 }
4427 }
4428
4429 MFloat blEdgeVValueGlobal = F0;
4430 MPI_Allreduce(&blEdgeVValueLocal, &blEdgeVValueGlobal, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_rescalingCommGrComm, AT_,
4431 "blEdgeVValueLocal", "blEdgeVValueGlobal");
4432
4433 if(edgePointIsSet) {
4434 for(MInt j = edgePointJ; j < jEnd; ++j) {
4435 const MInt cellId = cellIndex(i, j);
4436 m_cells->pvariables[PV->V][cellId] = blEdgeVValueGlobal;
4437 const MFloat pres = PV->PInfinity;
4438 m_cells->pvariables[PV->P][cellId] = pres;
4439 }
4440 }
4441
4442 for(MInt j = 0; j < m_nCells[0]; ++j) {
4443 // extrapolation for second GC
4444 const MInt cellId = cellIndex(1, j);
4445 const MInt cellIdM1 = cellIndex(0, j);
4446 const MInt cellIdadj = cellIndex(2, j);
4447
4448 for(MInt var = 0; var < PV->noVariables; var++) {
4449 m_cells->pvariables[var][cellIdM1] = 2.0 * m_cells->pvariables[var][cellId] - m_cells->pvariables[var][cellIdadj];
4450 }
4451 }
4452}
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
Solver::domainId
virtual MInt domainId() const
Return the domainId (rank)
Definition: solver.h:383
Solver::m_Ma
MFloat m_Ma
the Mach number
Definition: solver.h:71
mMin
constexpr T mMin(const T &x, const T &y)
Definition: functions.h:90
mMax
constexpr T mMax(const T &x, const T &y)
Definition: functions.h:94
MBool
bool MBool
Definition: maiatypes.h:58

◆ bc2511()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2511 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 4457 of file fvstructuredbndrycnd2d.cpp.

4457 {
4458 MInt* start = m_physicalBCMap[bcId]->start1;
4459 cout.precision(8);
4460
4461 // scaling between delta0 and delta2
4462 const MFloat F727 = 72.0 / 7.0; // 8.5, 8.0
4463 const MInt i = start[0];
4464 const MFloat yWall = F0; // this has been fixed else method does not works
4465 const MFloat maxIntegrationHeight = 2.0 * m_rescalingBLT;
4466
4470
4471 // thetaLocal.fill(F0); //initialize scratch space to zero // only for parallel use
4472 MFloatScratchSpace thetaLocal(2, AT_, "thetaLocalRe");
4473 MFloatScratchSpace thetaGlobal(2, AT_, "thetaGlobalRe");
4474 thetaLocal.fill(F0); // initialize scratch space
4475 thetaGlobal.fill(F0);
4476
4477 // the offest position in k-direction is the offset
4478 // compute the local moment thickness j=direction of integration
4479 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; ++j) {
4480 const MInt cellId = cellIndex(i, j);
4481 const MInt pointIdM1 = getPointIdFromCell(i, j);
4482 const MInt pointIdP1 = getPointIdFromPoint(pointIdM1, 0, 1);
4483
4484 if(m_grid->m_coordinates[1][pointIdM1] > maxIntegrationHeight) {
4485 continue;
4486 }
4487
4488 const MFloat urat = m_cells->pvariables[PV->U][cellId] / PV->UInfinity;
4489 const MFloat momThick = m_cells->pvariables[PV->U][cellId] * m_cells->pvariables[PV->RHO][cellId] * fabs(F1 - urat)
4490 / (CV->rhoUInfinity);
4491
4492 // integrate normal to the wall
4493 const MFloat ydist = m_grid->m_coordinates[1][pointIdP1] - m_grid->m_coordinates[1][pointIdM1];
4494 thetaLocal(1) += momThick * ydist;
4495 }
4496
4497
4498 // communicate the Thickness across the plane
4499 MPI_Allreduce(thetaLocal.begin(), thetaGlobal.begin(), 2, MPI_DOUBLE, MPI_SUM, m_solver->m_rescalingCommGrComm, AT_,
4500 "thetaLocal.begin()", "thetaGlobal.begin()");
4501
4505
4506 const MFloat delta = F727 * thetaGlobal(1);
4507 const MInt noVar = 3; // for more variables if wanted
4508 const MInt wallLocalOffset = m_solver->m_nOffsetCells[0]; // Offset in j-direction
4509
4510 MFloatScratchSpace wallPropertiesLocal(noVar, AT_, "wallPropertiesLocalRe");
4511 MFloatScratchSpace wallProperties(noVar, AT_, "wallPropertiesRe");
4512 wallPropertiesLocal.fill(F0);
4513 wallProperties.fill(F0);
4514
4515 // determine the wall stuff if wall is contained whithin the partition
4516 if(wallLocalOffset == 0 && m_solver->m_nActiveCells[1] >= m_noGhostLayers) {
4517 const MInt cellId = cellIndex(i, m_noGhostLayers);
4518 const MFloat rho = m_cells->pvariables[PV->RHO][cellId];
4519 const MFloat p = m_cells->pvariables[PV->P][cellId];
4520 const MFloat t = p * m_solver->m_gamma / rho;
4521 const MFloat mu = SUTHERLANDLAW(t);
4522 const MFloat uWall = fabs(m_cells->pvariables[PV->U][cellId]);
4523 const MFloat ydist = m_cells->coordinates[1][cellId] - yWall;
4524 const MFloat uTau = sqrt(uWall * mu / (ydist * rho));
4525
4526 wallPropertiesLocal(0) = uTau;
4527 wallPropertiesLocal(1) = rho;
4528 wallPropertiesLocal(2) = t;
4529 }
4530
4531 MPI_Allreduce(wallPropertiesLocal.begin(), wallProperties.begin(), noVar, MPI_DOUBLE, MPI_SUM,
4532 m_solver->m_rescalingCommGrComm, AT_, "wallPropertiesLocal.begin()", "wallProperties.begin()");
4533
4534 MFloat tempWallInletLocal = F0;
4535 MFloat tempWallInletGlobal = F0;
4536 MPI_Allreduce(&tempWallInletLocal, &tempWallInletGlobal, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_rescalingCommGrComm, AT_,
4537 "tempWallInletLocal", "tempWallInletGlobal");
4538
4542
4543 MInt totalCells = m_grid->getMyBlockNoCells(0) + 1;
4544
4545 const MInt noVariables = PV->noVariables + 1;
4546 MFloatScratchSpace varSliceLocal(noVariables, totalCells, AT_, "varSliceLocal");
4547 MFloatScratchSpace varSlice(noVariables, totalCells, AT_, "varSlice");
4548
4549 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; ++j) {
4550 const MInt cellId = cellIndex(i, j);
4551 const MFloat rho = m_cells->pvariables[PV->RHO][cellId];
4552 const MFloat frho = F1 / rho;
4553 const MFloat p = pressure(cellId);
4554 const MFloat temp = p * m_solver->m_gamma * frho;
4555 const MFloat mu = SUTHERLANDLAW(temp);
4556 const MFloat uTauRe = wallProperties(0);
4557 const MFloat yIn = (m_cells->coordinates[1][cellId] - yWall) * uTauRe * rho / (mu * sqrt(m_solver->m_Re0));
4558 const MFloat yOut = (m_cells->coordinates[1][cellId] - yWall) * rho / (delta * CV->rhoInfinity);
4559 const MFloat u = m_cells->pvariables[PV->U][cellId];
4560 const MFloat v = m_cells->pvariables[PV->V][cellId];
4561
4562 //>RANS
4563 const MFloat mut = m_cells->pvariables[PV->RANS_VAR[0]][cellId] * rho;
4564 //<RANS
4565
4566 const MInt localId = m_solver->m_nOffsetCells[0] + j - m_noGhostLayers + 1;
4567
4568 // save the variables u,v,w,t,yI,yO
4569 varSliceLocal(0, localId) = u;
4570 varSliceLocal(1, localId) = v;
4571 varSliceLocal(2, localId) = temp;
4572 varSliceLocal(3, localId) = yIn;
4573 varSliceLocal(4, localId) = yOut;
4574 varSliceLocal(5, localId) = mut;
4575 }
4576
4577 // set first value at the wall manually
4578 if(m_solver->m_nOffsetCells[0] == 0) {
4579 varSliceLocal(0, 0) = 0.0;
4580 varSliceLocal(1, 0) = 0.0;
4581 varSliceLocal(2, 0) = varSliceLocal(2, 1);
4582 varSliceLocal(3, 0) = 0.0;
4583 varSliceLocal(4, 0) = 0.0;
4584 varSliceLocal(5, 0) = varSliceLocal(5, 1);
4585 }
4586
4587 // communicate the slice
4588 MPI_Allreduce(varSliceLocal.begin(), varSlice.begin(), noVariables * totalCells, MPI_DOUBLE, MPI_SUM,
4589 m_solver->m_rescalingCommGrComm, AT_, "varSliceLocal.begin()", "varSlice.begin()");
4590
4591 MFloat blEdgeVValueLocal = F0;
4592 MFloat blEdgeVValueGlobal = F0;
4593 MPI_Allreduce(&blEdgeVValueLocal, &blEdgeVValueGlobal, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_rescalingCommGrComm, AT_,
4594 "blEdgeVValueLocal", "blEdgeVValueGlobal");
4595}

◆ bc2600()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2600 ( MInt  bcId)
virtual

Precribes a profile from the restart file extrapolate pressure from computational domain

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 4604 of file fvstructuredbndrycnd2d.cpp.

4604 {
4605 MInt* start = m_physicalBCMap[bcId]->start1;
4606 MInt* end = m_physicalBCMap[bcId]->end1;
4607
4608 switch(m_physicalBCMap[bcId]->face) {
4609 case 0: {
4610 for(MInt i = start[0]; i < end[0]; i++) {
4611 for(MInt j = start[1]; j < end[1]; j++) {
4612 MInt cellId = cellIndex(m_noGhostLayers - 1 - i, j);
4613 MInt cellIdadj = cellIndex(m_noGhostLayers - i, j);
4614
4615 // extrapolate pressure to ghost cells
4616 m_cells->pvariables[PV->P][cellId] = m_cells->pvariables[PV->P][cellIdadj];
4617 }
4618 }
4619 break;
4620 }
4621 default: {
4622 mTerm(1, AT_, "Face not implemented");
4623 break;
4624 }
4625 }
4626}

◆ bc2999()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc2999 ( MInt  bcId)
virtual

Prescribes the velocity distribution for a Blasius boundary layer at the inflow, with the Blasius helper functions

Author
Marian Albers

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2623 of file fvstructuredbndrycnd2d.cpp.

2623 {
2624 const MInt IJ[nDim] = {1, m_nCells[1]};
2625 MInt* start = m_physicalBCMap[bcId]->start1;
2626 MInt* end = m_physicalBCMap[bcId]->end1;
2627
2628 // Here we find out the normal direction of the
2629 // boundary and the two tangential directions.
2630 // This way we can make a general formulation of
2631 // the boundary condition
2632 const MInt face = m_physicalBCMap[bcId]->face;
2633 const MInt normalDir = face / 2;
2634 const MInt firstTangentialDir = (normalDir + 1) % nDim;
2635 const MInt normalDirStart = start[normalDir];
2636 const MInt firstTangentialStart = start[firstTangentialDir];
2637 const MInt firstTangentialEnd = end[firstTangentialDir];
2638 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
2639 // determine indices for direction help
2640 const MInt n = (face % 2) * 2 - 1; //-1,+1
2641 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
2642 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
2643 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
2644 // const MInt a2 = normalDirStart + (MInt)(0.5-(2.5*(MFloat)n)); //+3,-2
2645
2646 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
2647 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
2648 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
2649 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
2650 // const MInt cellIdA2 = a2*inc[0] + t1*inc[1] + t2*inc[2];
2651
2652 MFloat vel[nDim];
2653 m_solver->getBlasiusVelocity(cellIdG1, vel);
2654 m_cells->pvariables[PV->U][cellIdG1] = vel[0];
2655 m_cells->pvariables[PV->V][cellIdG1] = vel[1];
2656 m_cells->pvariables[PV->P][cellIdG1] = m_cells->pvariables[PV->P][cellIdA1];
2657 m_cells->pvariables[PV->RHO][cellIdG1] = CV->rhoInfinity;
2658
2659 m_solver->getBlasiusVelocity(cellIdG2, vel);
2660 m_cells->pvariables[PV->U][cellIdG2] = vel[0];
2661 m_cells->pvariables[PV->V][cellIdG2] = vel[1];
2662 m_cells->pvariables[PV->P][cellIdG2] =
2663 F2 * m_cells->pvariables[PV->P][cellIdG1] - m_cells->pvariables[PV->P][cellIdA1];
2664 m_cells->pvariables[PV->RHO][cellIdG2] = CV->rhoInfinity;
2665 }
2666}
FvStructuredSolver::getBlasiusVelocity
virtual void getBlasiusVelocity(MInt cellId, MFloat *const vel)
Load variables for the specified timeStep.
Definition: fvstructuredsolver.cpp:7721

◆ bc3000()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc3000 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 3686 of file fvstructuredbndrycnd2d.cpp.

3686 {
3687 TRACE();
3688 MInt* start = m_physicalBCMap[bcId]->start1;
3689 MInt* end = m_physicalBCMap[bcId]->end1;
3690
3691 switch(m_physicalBCMap[bcId]->face) {
3692 case 0: {
3693 mTerm(1, "Not implemted yet!");
3694 break;
3695 }
3696 case 1: {
3697 mTerm(1, "Not implemted yet!");
3698 break;
3699 }
3700 case 2: {
3701 MInt cellId = -1;
3702 MInt cellIdadj = -1;
3703 const MInt cellShift = 2 * m_noGhostLayers - 1;
3704
3705 for(MInt j = start[1]; j < end[1]; j++) {
3706 for(MInt i = start[0]; i < end[0]; i++) {
3707 cellId = cellIndex(i, j); // ghost
3708 cellIdadj = cellIndex(i, cellShift - j); // field
3709
3710 m_cells->pvariables[PV->RHO][cellId] = m_cells->pvariables[PV->RHO][cellIdadj];
3711 m_cells->pvariables[PV->P][cellId] = m_cells->pvariables[PV->P][cellIdadj];
3712 m_cells->pvariables[PV->U][cellId] = m_cells->pvariables[PV->U][cellIdadj];
3713 m_cells->pvariables[PV->V][cellId] = -m_cells->pvariables[PV->V][cellIdadj];
3714 if(isRans) {
3715 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3716 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] =
3717 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdadj];
3718 }
3719 }
3720 }
3721 }
3722
3723 break;
3724 }
3725 case 3: {
3726 MInt cellId = -1;
3727 MInt cellIdadj = -1;
3728 const MInt cellShift = 2 * (m_nCells[0] - 1) - 2 * m_noGhostLayers + 1;
3729
3730 for(MInt j = start[1]; j < end[1]; j++) {
3731 for(MInt i = start[0]; i < end[0]; i++) {
3732 // mirroring
3733 cellId = cellIndex(i, j); // ghost
3734 cellIdadj = cellIndex(i, cellShift - j); // field
3735
3736 m_cells->pvariables[PV->RHO][cellId] = m_cells->pvariables[PV->RHO][cellIdadj];
3737 m_cells->pvariables[PV->P][cellId] = m_cells->pvariables[PV->P][cellIdadj];
3738 m_cells->pvariables[PV->U][cellId] = m_cells->pvariables[PV->U][cellIdadj];
3739 m_cells->pvariables[PV->V][cellId] = -m_cells->pvariables[PV->V][cellIdadj];
3740 if(isRans) {
3741 for(MInt ransVarId = 0; ransVarId < m_solver->m_noRansEquations; ++ransVarId) {
3742 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellId] =
3743 m_cells->pvariables[PV->RANS_VAR[ransVarId]][cellIdadj];
3744 }
3745 }
3746 }
3747 }
3748 break;
3749 }
3750
3751 default: {
3752 mTerm(1, AT_, "Face direction not implemented)");
3753 }
3754 }
3755}

◆ bc6002()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::bc6002 ( MInt  bcId)
overridevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 4632 of file fvstructuredbndrycnd2d.cpp.

4632 {
4633 TRACE();
4634
4635 const MInt IJ[nDim] = {1, m_nCells[1]};
4636 MInt* start = m_physicalBCMap[bcId]->start1;
4637 MInt* end = m_physicalBCMap[bcId]->end1;
4638
4639 constexpr MFloat eps = 1e-8;
4640 const MFloat gamma = m_solver->m_gamma;
4641 const MFloat gammaOverGammaMinusOne = gamma / (gamma - 1);
4642
4643 if(m_physicalBCMap[bcId]->Nstar == -1) {
4644 // Here we find out the normal direction of the
4645 // boundary and the tangential direction.
4646 // This way we can make a general formulation of
4647 // the boundary condition
4648 const MInt face = m_physicalBCMap[bcId]->face;
4649 const MInt normalDir = face / 2;
4650 const MInt firstTangentialDir = (normalDir + 1) % nDim;
4651 const MInt normalDirStart = start[normalDir];
4652 const MInt firstTangentialStart = start[firstTangentialDir];
4653 const MInt firstTangentialEnd = end[firstTangentialDir];
4654 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
4655 // determine indices for direction help
4656 const MInt n = (face % 2) * 2 - 1; //-1,+1
4657 const MInt g[2] = {normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)), //+1,0
4658 normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n))}; // 0,+1
4659 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
4660
4661 // convention: index x is unknown and y is known
4662 for(MInt t1 = firstTangentialStart, pos = 0; t1 < firstTangentialEnd; t1++, ++pos) {
4663 const MInt cellIdG1 = g[0] * inc[0] + t1 * inc[1];
4664 const MFloat normalVec[nDim] = {m_cells->fq[FQ->NORMAL[0]][cellIdG1], m_cells->fq[FQ->NORMAL[1]][cellIdG1]};
4665 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
4666 const MFloat por_x = m_cells->fq[FQ->POROSITY][cellIdA1];
4667 for(MInt i = 0; i < 2 /*m_noGhostLayers*/; ++i) {
4668 const MInt cellIdG = g[i] * inc[0] + t1 * inc[1];
4669
4670 const MFloat por_y = m_cells->fq[FQ->POROSITY][cellIdG];
4671 const MFloat rho_y = m_cells->pvariables[PV->RHO][cellIdG];
4672 const MFloat p_y = m_cells->pvariables[PV->P][cellIdG];
4673 const MFloat u_y = m_cells->pvariables[PV->U][cellIdG];
4674 const MFloat v_y = m_cells->pvariables[PV->V][cellIdG];
4675 const MFloat k_y = (isRans) ? m_cells->pvariables[PV->RANS_VAR[0]][cellIdG] : 0.0;
4676
4677 // Compute normal and tangential velocity components
4678 const MFloat U_y = u_y * normalVec[0] + v_y * normalVec[1];
4679 const MFloat V_y[nDim] = {u_y - U_y * normalVec[0], v_y - U_y * normalVec[1]};
4680
4681 //
4682 const MFloat auxconst1 = gammaOverGammaMinusOne * p_y / pow(rho_y, gamma);
4683 const MFloat auxconst2 = (por_y / por_x) * rho_y * k_y;
4684 const MFloat auxconst3 = -(gammaOverGammaMinusOne * p_y / rho_y + 0.5 * POW2(U_y) + k_y);
4685 const MFloat auxconst4 = 0.5 * POW2(por_y / por_x * rho_y * U_y);
4686
4687 // TODO_SS labels:FV Later take a more sufficient initial value
4688 MFloat rho_x = rho_y;
4689
4690 MFloat res = auxconst1 * pow(rho_x, gamma + 1) + auxconst2 * rho_x + auxconst3 * POW2(rho_x) + auxconst4;
4691
4692 MInt testcounter = 0;
4693 while(res > eps) {
4694 // TODO_SS labels:FV,totest Check if it is divided by zero
4695 const MFloat dresdrho = (gamma + 1) * auxconst1 * pow(rho_x, gamma) + auxconst2 + 2 * auxconst3 * rho_x;
4696 const MFloat deltarho_x = -res / dresdrho;
4697 rho_x += deltarho_x;
4698
4699 res = auxconst1 * pow(rho_x, gamma + 1) + auxconst2 * rho_x + auxconst3 * POW2(rho_x) + auxconst4;
4700
4701 ++testcounter;
4702 if(testcounter > 100) break;
4703 }
4704
4705 if(testcounter > 10) {
4706 cout << "Loop in BC6002 took " << testcounter << " steps!"
4707 << " domainId=" << m_solver->domainId() << " cellId: " << cellIdG << " (" << g[i] << "|" << t1
4708 << ") res=" << res << " (x|y)=" << m_cells->coordinates[0][cellIdG] << "|"
4709 << m_cells->coordinates[1][cellIdG] << por_y << " " << rho_y << " " << p_y << " " << u_y << " " << v_y
4710 << endl;
4711 }
4712
4713 // Compute remaining variables
4714 const MFloat temp = por_y * rho_y / (por_x * rho_x);
4715 const MFloat p_x = p_y * pow(rho_x / rho_y, gamma);
4716 const MFloat U_x = temp * U_y;
4717 const MFloat u_x = U_x * normalVec[0] + V_y[0] /*=V_x*/;
4718 const MFloat v_x = U_x * normalVec[1] + V_y[1] /*=V_x*/;
4719
4720 // Assign solution to pvariables
4721 m_cells->pvariables[PV->RHO][cellIdG] = rho_x;
4722 m_cells->pvariables[PV->P][cellIdG] = p_x;
4723 m_cells->pvariables[PV->U][cellIdG] = u_x;
4724 m_cells->pvariables[PV->V][cellIdG] = v_x;
4725 if(isRans) {
4726 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG] *= temp;
4727 m_cells->pvariables[PV->RANS_VAR[1]][cellIdG] *= temp;
4728 }
4729 }
4730 }
4731 } else {
4732 const MInt singularId = m_physicalBCMap[bcId]->SingularId;
4733 const auto& singularity = m_solver->m_singularity[singularId];
4734
4735 // Check
4736 const MInt* start2 = singularity.start;
4737 const MInt* end2 = singularity.end;
4738 for(MInt d = 0; d < nDim; ++d) {
4739 if(end[d] - start[d] != end2[d] - start2[d]) mTerm(1, "SCHEISSE");
4740 }
4741
4742 // cout << globalTimeStep<< ": SINGULARITY BC d=" << m_solver->domainId() << " start=" << start[0] << "|" <<
4743 // start[1] << " start2=" << start2[0] << "|" << start2[1] << " end=" << end[0] << "|" << end[1] << " end2=" <<
4744 // end2[0] << "|" << end2[1] << endl;
4745
4746 for(MInt j = start[1], jj = start2[1]; j < end[1]; j++, jj++) {
4747 for(MInt i = start[0], ii = start2[0]; i < end[0]; i++, ii++) {
4748 const MInt cellIdG = cellIndex(i, j);
4749 const MInt cellIdA1 = cellIndex(ii, jj);
4750 const MFloat normalVec[nDim] = {m_cells->fq[FQ->NORMAL[0]][cellIdG], m_cells->fq[FQ->NORMAL[1]][cellIdG]};
4751 const MFloat por_x = m_cells->fq[FQ->POROSITY][cellIdA1];
4752
4753 const MFloat por_y = m_cells->fq[FQ->POROSITY][cellIdG];
4754 const MFloat rho_y = m_cells->pvariables[PV->RHO][cellIdG];
4755 const MFloat p_y = m_cells->pvariables[PV->P][cellIdG];
4756 const MFloat u_y = m_cells->pvariables[PV->U][cellIdG];
4757 const MFloat v_y = m_cells->pvariables[PV->V][cellIdG];
4758 const MFloat k_y = (isRans) ? m_cells->pvariables[PV->RANS_VAR[0]][cellIdG] : 0.0;
4759
4760 // Compute normal and tangential velocity components
4761 const MFloat U_y = u_y * normalVec[0] + v_y * normalVec[1];
4762 const MFloat V_y[nDim] = {u_y - U_y * normalVec[0], v_y - U_y * normalVec[1]};
4763
4764 //
4765 const MFloat auxconst1 = gammaOverGammaMinusOne * p_y / pow(rho_y, gamma);
4766 const MFloat auxconst2 = (por_y / por_x) * rho_y * k_y;
4767 const MFloat auxconst3 = -(gammaOverGammaMinusOne * p_y / rho_y + 0.5 * POW2(U_y) + k_y);
4768 const MFloat auxconst4 = 0.5 * POW2(por_y / por_x * rho_y * U_y);
4769
4770 // TODO_SS labels:FV Later take a more sufficient initial value
4771 MFloat rho_x = rho_y;
4772
4773 MFloat res = auxconst1 * pow(rho_x, gamma + 1) + auxconst2 * rho_x + auxconst3 * POW2(rho_x) + auxconst4;
4774
4775 MInt testcounter = 0;
4776 while(res > eps) {
4777 // TODO_SS labels:FV,totest Check if it is divided by zero
4778 const MFloat dresdrho = (gamma + 1) * auxconst1 * pow(rho_x, gamma) + auxconst2 + 2 * auxconst3 * rho_x;
4779 const MFloat deltarho_x = -res / dresdrho;
4780 rho_x += deltarho_x;
4781
4782 res = auxconst1 * pow(rho_x, gamma + 1) + auxconst2 * rho_x + auxconst3 * POW2(rho_x) + auxconst4;
4783
4784 ++testcounter;
4785 if(testcounter > 100) break;
4786 }
4787
4788 if(testcounter > 10) {
4789 cout << "Loop in BC6002 took " << testcounter << " steps!"
4790 << " domainId=" << m_solver->domainId() << " cellId: " << cellIdG << " res=" << res
4791 << " (x|y)=" << m_cells->coordinates[0][cellIdG] << "|" << m_cells->coordinates[1][cellIdG] << por_y
4792 << " " << rho_y << " " << p_y << " " << u_y << " " << v_y << endl;
4793 }
4794
4795 // Compute remaining variables
4796 const MFloat temp = por_y * rho_y / (por_x * rho_x);
4797 const MFloat p_x = p_y * pow(rho_x / rho_y, gamma);
4798 const MFloat U_x = temp * U_y;
4799 const MFloat u_x = U_x * normalVec[0] + V_y[0] /*=V_x*/;
4800 const MFloat v_x = U_x * normalVec[1] + V_y[1] /*=V_x*/;
4801
4802 // Assign solution to pvariables
4803 m_cells->pvariables[PV->RHO][cellIdG] = rho_x;
4804 m_cells->pvariables[PV->P][cellIdG] = p_x;
4805 m_cells->pvariables[PV->U][cellIdG] = u_x;
4806 m_cells->pvariables[PV->V][cellIdG] = v_x;
4807 if(isRans) {
4808 m_cells->pvariables[PV->RANS_VAR[0]][cellIdG] *= temp;
4809 m_cells->pvariables[PV->RANS_VAR[1]][cellIdG] *= temp;
4810 }
4811 }
4812 }
4813 }
4814}
FvStructuredSolver::m_singularity
SingularInformation * m_singularity
Definition: fvstructuredsolver.h:835
StructuredBndryCnd< 2 >::FQ
std::unique_ptr< StructuredFQVariables > & FQ
Definition: fvstructuredbndrycnd.h:55

◆ calc_cp_cf()

template<MBool isRans>
template<MBool calcCp, MBool calcCf, MBool interface>
template void StructuredBndryCnd2D< isRans >::calc_cp_cf< true, true, true > ( const  MInt,
const  MInt,
const  MInt,
const  MInt,
MFloat(&)  [calcCp+nDim *calcCf] 
)

Definition at line 2133 of file fvstructuredbndrycnd2d.cpp.

2134 {
2135 // cf=nu*du/dy/(rho*u_8*u_8*0.5)
2136 // cp=2*(p-p_8)/(rho8*u_8*u_8*0.5)
2137
2138 static const MFloat UT = m_solver->m_Ma * sqrt(PV->TInfinity);
2139 static const MFloat fstagnationPressure = F1 / (CV->rhoInfinity * F1B2 * POW2(UT));
2140 static const MFloat fre0 = F1 / (m_solver->m_Re0);
2141
2142 // first get the position of the wall (reference!!)
2143 // x is always used as coordinate in normal direction
2144 const MFloat xRef = 0.5 * (m_grid->m_coordinates[0][pIJ] + m_grid->m_coordinates[0][pIJP]);
2145 const MFloat yRef = 0.5 * (m_grid->m_coordinates[1][pIJ] + m_grid->m_coordinates[1][pIJP]);
2146
2147 // compute the pressure coefficient
2148 //-->
2149 // get the distcance
2150 const MFloat dx2 = sqrt(POW2(m_cells->coordinates[0][cellIdP1] - m_cells->coordinates[0][cellId])
2151 + POW2(m_cells->coordinates[1][cellIdP1] - m_cells->coordinates[1][cellId]));
2152 const MFloat dx1 = sqrt(POW2(m_cells->coordinates[0][cellId] - xRef) + POW2(m_cells->coordinates[1][cellId] - yRef));
2153
2154 if(calcCp) {
2155 // pressures
2156 const MFloat p1 = m_cells->pvariables[PV->P][cellId];
2157 const MFloat p2 = m_cells->pvariables[PV->P][cellIdP1];
2158 // extrapolation to the face
2159 const MFloat pW = ((p1 - p2) / dx2) * dx1 + p1;
2160 const MFloat cpn = (pW - PV->PInfinity) * fstagnationPressure;
2161 output[0] = cpn;
2162 //<--- finished computation of cp
2163 }
2164
2165 if(calcCf) {
2166 // compute the skin-friction coefficient
2167 //-->
2168 // compute orthogonal distance surface-plane to cell center
2169 MFloat supportVec[nDim] = {};
2170 MFloat firstVec[nDim] = {};
2171 MFloat normalVec[nDim] = {};
2172 MFloat cellVec[nDim] = {};
2173 for(MInt dim = 0; dim < nDim; dim++) {
2174 supportVec[dim] = m_grid->m_coordinates[dim][pIJ];
2175 firstVec[dim] = m_grid->m_coordinates[dim][pIJP] - m_grid->m_coordinates[dim][pIJ];
2176 cellVec[dim] = m_cells->coordinates[dim][cellId];
2177 }
2178
2179 normalVec[0] = -firstVec[1];
2180 normalVec[1] = firstVec[0];
2181 const MFloat normalLength = sqrt(POW2(normalVec[0]) + POW2(normalVec[1]));
2182 MFloat orthDist = F0;
2183
2184 for(MInt dim = 0; dim < nDim; dim++) {
2185 orthDist += (cellVec[dim] - supportVec[dim]) * normalVec[dim];
2186 }
2187
2188 orthDist = fabs(orthDist / normalLength);
2189
2190 // extrapolate the Temperature to the wall
2191 const MFloat T1 = temperature(cellId);
2192 const MFloat T2 = temperature(cellIdP1);
2193
2194 // extrapolation to the face
2195 const MFloat tBc = ((T1 - T2) / dx2) * dx1 + T1;
2196 const MFloat mue = SUTHERLANDLAW(tBc);
2197
2198 // velocities
2199 const MFloat u1 = m_cells->pvariables[PV->U][cellId];
2200 const MFloat v1 = m_cells->pvariables[PV->V][cellId];
2201
2202 // at 6000er-interfaces the velocity is not zero
2203 MFloat uBc = 0.0;
2204 MFloat vBc = 0.0;
2205 if(interface) {
2206 const MFloat u2 = m_cells->pvariables[PV->U][cellIdP1];
2207 const MFloat v2 = m_cells->pvariables[PV->V][cellIdP1];
2208 uBc = ((u1 - u2) / dx2) * dx1 + u1;
2209 vBc = ((v1 - v2) / dx2) * dx1 + v1;
2210 } else {
2211 if(m_solver->m_movingGrid) {
2212 uBc = F1B2 * (m_grid->m_velocity[0][pIJ] + m_grid->m_velocity[0][pIJP]);
2213 vBc = F1B2 * (m_grid->m_velocity[1][pIJ] + m_grid->m_velocity[1][pIJP]);
2214 }
2215 }
2216
2217 // compute gradient
2218 MFloat dudeta = (u1 - uBc) / orthDist;
2219 MFloat dvdeta = (v1 - vBc) / orthDist;
2220
2221 // using taylor expansion a second-order approximation
2222 // of du/dn can be achieved at the wall
2223 if(m_solver->m_forceSecondOrder) {
2224 const MFloat u2 = m_cells->pvariables[PV->U][cellIdP1];
2225 const MFloat v2 = m_cells->pvariables[PV->V][cellIdP1];
2226 const MFloat dx3 = dx1 + dx2;
2227 dudeta = (u1 * dx3 * dx3 + (dx1 * dx1 - dx3 * dx3) * uBc - u2 * dx1 * dx1) / (dx1 * dx3 * dx3 - dx1 * dx1 * dx3);
2228 dvdeta = (v1 * dx3 * dx3 + (dx1 * dx1 - dx3 * dx3) * vBc - v2 * dx1 * dx1) / (dx1 * dx3 * dx3 - dx1 * dx1 * dx3);
2229 }
2230
2231 const MFloat taux = mue * dudeta * fre0;
2232 const MFloat tauy = mue * dvdeta * fre0;
2233
2234 const MFloat cfx = taux * fstagnationPressure;
2235 const MFloat cfy = tauy * fstagnationPressure;
2236
2237
2238 output[calcCp + 0] = cfx, output[calcCp + 1] = cfy;
2239 }
2240}

◆ cellIndex()

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::cellIndex ( MInt  i,
MInt  j 
)
inline

Definition at line 2500 of file fvstructuredbndrycnd2d.cpp.

2500 {
2501 return i + (j * m_nCells[1]);
2502}

◆ computeDistance2Map()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::computeDistance2Map ( const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &  wallDistInfo,
MFloat * const  r_cells_fq_distance,
std::vector< std::pair< MInt, MInt > > &  r_cellId2globalWallId,
std::vector< MInt > &  r_displs,
std::vector< MFloat > &  r_weighting 
)

Computes shortest distance to wall/fluid-porous interface, which is described by the grid as the sum of linear line elements and stores the information required to exchange interpolated data from the nearest neighbors to a cell.

Date
2020-03-09
Parameters
[in]wallDistInfoe.g maps of all walls or fluid-porous-interfaces
[out]r_cells_fq_distancee.g. m_cells->fq[FQ->WALLDISTANCE]
[out]r_cellId2globalWallId\parma[out] r_displs
Date
2020-03-09
Parameters
[in]wallDistInfoe.g maps of all walls or fluid-porous-interfaces
[out]r_cells_fq_distancee.g. m_cells->fq[FQ->WALLDISTANCE]
[out]r_cellId2globalWallId\parma[out] r_displs

Set distance, nearest neigbor cell identifiers & weighting

Definition at line 1366 of file fvstructuredbndrycnd2d.cpp.

1371 {
1372 // Initialize array with high numbers
1373 MFloat maxWallDistance = numeric_limits<MFloat>::max();
1374 maxWallDistance = Context::getSolverProperty<MFloat>("maxWallDistance", m_solverId, AT_, &maxWallDistance);
1375 for(MInt id = 0; id < m_noCells; ++id) {
1376 r_cells_fq_distance[id] = maxWallDistance;
1377 }
1378
1379 // const vector<StructuredWindowMap*> wallDistInfo = m_auxDataMap;//m_solver->m_windowInfo->m_wallDistInfoMap;
1380 const MInt noWallDistInfo = wallDistInfo.size(); // contains the size of the maps
1381 std::vector<MInt> wallNormalIndex(noWallDistInfo);
1382 std::vector<MInt> normalDir(noWallDistInfo);
1383 std::vector<std::array<MInt, 2 * nDim>> startEndTangential(noWallDistInfo);
1384 // Number of grid points on current rank to store
1385 MInt cellmemSize = 0;
1386 for(MInt nn = 0; nn < noWallDistInfo; ++nn) {
1387 // start & end represent exactly the range of grid point indices, which form the map/boundary
1388 const MInt* const start = wallDistInfo[nn]->start1;
1389 const MInt* const end = wallDistInfo[nn]->end1;
1390 const MInt face = wallDistInfo[nn]->face;
1391 normalDir[nn] = face / 2;
1392 wallNormalIndex[nn] = start[normalDir[nn]];
1393
1394 MInt size = 1;
1395 for(MInt dim = 0; dim < nDim - 1; ++dim) {
1396 const MInt tangentialDir = (normalDir[nn] + (dim + 1)) % nDim;
1397 startEndTangential[nn][dim * 2 + 0] = start[tangentialDir]; // + m_noGhostLayers;
1398 startEndTangential[nn][dim * 2 + 1] =
1399 end[tangentialDir] + 1; //+1, because grid coords not cell coords// - m_noGhostLayers;
1400 size *= startEndTangential[nn][dim * 2 + 1] - startEndTangential[nn][dim * 2 + 0];
1401 }
1402 cellmemSize += size;
1403#ifndef NDEBUG
1404 cout << "Debug output: wallDistInfo=" << nn << ": face=" << face << "; normalDir=" << normalDir[nn]
1405 << "; wallNormalIndex=" << wallNormalIndex[nn] << "; startEndTangential=" << startEndTangential[nn][0] << "|"
1406 << startEndTangential[nn][1] << "normal start/end=" << start[normalDir[nn]] << "|" << end[normalDir[nn]]
1407 << endl;
1408#endif
1409 }
1410
1411 // 2.2) allocate the space for all the coordinates in Scratch!
1412 // memory will not be needed later!
1413 MFloatScratchSpace coordGridMem(std::max(1, nDim * cellmemSize), AT_, "coordGridMem");
1414 MFloatPointerScratchSpace wallDistSurfCoords(nDim, AT_, "wallDistSurfCoords");
1415 for(MInt dim = 0; dim < nDim; ++dim)
1416 wallDistSurfCoords[dim] = &coordGridMem[dim * cellmemSize];
1417 // Store if a grid point has a neighbor to the left or right on the same domain
1418 MIntScratchSpace isStartEndPoint(std::max(1, cellmemSize), AT_, "isStartEndPoint");
1419
1420 // 4) Put wall grid points and information about existence of neighbor grid point in temporary buffers;
1421 MInt cnt = 0;
1422 for(MInt nn = 0; nn < noWallDistInfo; ++nn) {
1423 MInt tIdx;
1424 MInt& i = (normalDir[nn] == 0) ? wallNormalIndex[nn] : tIdx;
1425 MInt& j = (normalDir[nn] == 0) ? tIdx : wallNormalIndex[nn];
1426 for(tIdx = startEndTangential[nn][0]; tIdx < startEndTangential[nn][1]; ++tIdx) {
1427 const MInt p = getPointIdFromCell(i, j);
1428 wallDistSurfCoords[0][cnt] = m_grid->m_coordinates[0][p];
1429 wallDistSurfCoords[1][cnt] = m_grid->m_coordinates[1][p];
1430 // WRONG: The case no neighbor to the left and right can not occur
1431 if(tIdx == startEndTangential[nn][0] && tIdx == startEndTangential[nn][1] - 1)
1432 isStartEndPoint[cnt] = 3;
1433 else if(tIdx == startEndTangential[nn][0])
1434 isStartEndPoint[cnt] = 1; // no neighbor to the left on this rank
1435 else if(tIdx == startEndTangential[nn][1] - 1)
1436 isStartEndPoint[cnt] = -1; // no neighbor to the right on this rank
1437 else
1438 isStartEndPoint[cnt] =
1439 0; // TODO_SS labels:FV,GRID,totest is this necessary, or will it be default initialized anyway?
1440 ++cnt;
1441 }
1442 }
1443
1444 // Broadcast all wall grid coordinates and isStartEndPoint to all ranks
1445 MInt noRanks;
1446 MPI_Comm_size(m_solver->m_StructuredComm, &noRanks);
1447 // recvcounts: How many wall grid coordinates to receive from each domain
1448 std::vector<MInt> recvcounts(noRanks);
1449 MPI_Allgather(&cellmemSize, 1, MPI_INT, recvcounts.data(), 1, MPI_INT, m_solver->m_StructuredComm, AT_, "cellmemSize",
1450 "recvcounts");
1451 r_displs.resize(noRanks);
1452 for(MInt r = 0; r < noRanks - 1; ++r) {
1453 r_displs[r + 1] = r_displs[r] + recvcounts[r];
1454 }
1455 const MInt noTotalWallPoints = r_displs[noRanks - 1] + recvcounts[noRanks - 1];
1456 // wallDistSurfCoordsAll: contains the wall grid coordinates of all ranks
1457 std::vector<std::vector<MFloat>> wallDistSurfCoordsAll(nDim);
1458 for(MInt dim = 0; dim < nDim; ++dim) {
1459 wallDistSurfCoordsAll[dim].resize(noTotalWallPoints);
1460 MPI_Allgatherv(&wallDistSurfCoords[dim][0], cellmemSize, MPI_DOUBLE, wallDistSurfCoordsAll[dim].data(),
1461 recvcounts.data(), r_displs.data(), MPI_DOUBLE, m_solver->m_StructuredComm, AT_,
1462 "wallDistSurfCoords", "wallDistSurfCoordsAll");
1463 }
1464 std::vector<MInt> isStartEndPointAll(noTotalWallPoints);
1465 MPI_Allgatherv(&isStartEndPoint[0], cellmemSize, MPI_INT, isStartEndPointAll.data(), recvcounts.data(),
1466 r_displs.data(), MPI_INT, m_solver->m_StructuredComm, AT_, "isStartEndPoint", "isStartEndPointAll");
1467 // Also send own inputBoxId
1468 std::vector<MInt> blockIds(noRanks, -1);
1469 const MInt myBlockId = m_solver->m_blockId;
1470 MPI_Allgather(&myBlockId, 1, MPI_INT, blockIds.data(), 1, MPI_INT, m_solver->m_StructuredComm, AT_, "myBlockId",
1471 "blockId");
1472
1473 // --> Now all ranks have all wall grid coordinates and the information if that grid point has any
1474 // neighbor to the left or right on same domain
1475
1476 // r_cellId2globalWallId: mapping of cellId to the positions of the two cells which encloses the nearest point at the
1477 // wall
1478 r_cellId2globalWallId.assign(m_noCells, std::make_pair(-1, -1));
1479 // weighting: weighting of the wall cell which is closest to a given cell; the weighting of its neighbor is
1480 // 1-weighting
1481 r_weighting.assign(m_noCells, 0.0);
1482
1483 // Return early: the return is done here because we need to allocate the vector
1484 // r_cellId2globalWallId, before we can exit this function
1485 if(noTotalWallPoints == 0) return;
1486
1487 // Build k-d-tree for a quick search
1488 // In the following Point<3> is used, but 3rd dimension is just dummy
1489 vector<Point<3>> pts;
1490 // cellIdShifts is used to retrieve the correct cellIds from the grid ids
1491 std::vector<MInt> cellIdShifts(noTotalWallPoints + 1, 0);
1492 for(MInt globalWallId = 0; globalWallId < noTotalWallPoints; ++globalWallId) {
1493 Point<3> a(wallDistSurfCoordsAll[0][globalWallId], wallDistSurfCoordsAll[1][globalWallId], 0, globalWallId);
1494 pts.push_back(a);
1495 cellIdShifts[globalWallId + 1] = cellIdShifts[globalWallId] - 1 * (isStartEndPointAll[globalWallId] == -1);
1496 }
1497 // build up the tree
1498 KDtree<3> tree(pts);
1499
1500 auto getDomainId = [&r_displs, noRanks](const MInt id, const MInt hint = 0) {
1501 for(MInt rank = hint; rank < (signed)r_displs.size() - 1; ++rank) {
1502 if(r_displs[rank + 1] > id) return rank;
1503 }
1504 return noRanks - 1;
1505 };
1506
1507 static constexpr MFloat radius = std::numeric_limits<MFloat>::max();
1508 static constexpr MFloat eps = 1e-16; // I hope this is a suitable value
1509
1510 // go through all the cells an determine the closest distance
1511 MInt results[3];
1512 MFloat dists[3];
1513 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1514 Point<3> pt(m_cells->coordinates[0][cellId], m_cells->coordinates[1][cellId], 0);
1515 MBool overflow = false;
1516 MInt nfound = tree.locatenearest(pt, radius, results, dists, 3, overflow);
1517 if(nfound != 3) mTerm(1, "Come on, your mesh should have at least 3 wall grid points!");
1518 MInt globalWallId1 = results[0];
1519 MInt globalWallId1_ = results[1];
1520 MInt globalWallId1__ = results[2];
1521
1522
1523 nfound = 1;
1524 const MFloat r2 = POW2(wallDistSurfCoordsAll[0][globalWallId1] - wallDistSurfCoordsAll[0][globalWallId1_])
1525 + POW2(wallDistSurfCoordsAll[1][globalWallId1] - wallDistSurfCoordsAll[1][globalWallId1_]);
1526 const MFloat r3 = POW2(wallDistSurfCoordsAll[0][globalWallId1] - wallDistSurfCoordsAll[0][globalWallId1__])
1527 + POW2(wallDistSurfCoordsAll[1][globalWallId1] - wallDistSurfCoordsAll[1][globalWallId1__]);
1528 nfound += r2 < eps;
1529 nfound += r3 < eps;
1530 // It might happen that all three points have same distance to target point, but point 1 and 3 coincide and not
1531 // 1 and 2
1532 if(nfound == 2 && r3 < r2) {
1533 const MInt temp = globalWallId1_;
1534 globalWallId1_ = globalWallId1__;
1535 globalWallId1__ = temp;
1536 }
1537 // Very special case: Because of the domain decomposition the start or end of a wall is on one domain without a
1538 // neighbor
1539 if(nfound == 2) {
1540 if(isStartEndPointAll[globalWallId1] == 3 && isStartEndPointAll[globalWallId1_] == 3) mTerm(1, "");
1541 if(isStartEndPointAll[globalWallId1] == 3) {
1542 globalWallId1 = globalWallId1_;
1543 nfound = 1;
1544 }
1545 if(isStartEndPointAll[globalWallId1_] == 3) {
1546 nfound = 1;
1547 }
1548 }
1549
1553 // Case 1: wall grid point has only one neighbor, because it is the begin or end of a solid wall
1554 // Case 2: wall grid point has two neighbors, both on one rank (nfound==1)
1555 // Case 3: wall grid point has two neighbors, but one neighbor is on a different rank (nfound==2)
1556 MInt globalCellId1 = -1;
1557 MInt globalCellId2 = -1;
1558 MFloat distance = std::numeric_limits<MFloat>::max(), weight;
1559 if(nfound == 1 && abs(isStartEndPointAll[globalWallId1]) == 1) {
1560 // Case 1:
1561 const MInt globalWallId2 = globalWallId1 + isStartEndPointAll[globalWallId1];
1562 globalCellId1 = (isStartEndPointAll[globalWallId1] == -1) ? globalWallId1 - 1 : globalWallId1;
1563 globalCellId1 += cellIdShifts[globalWallId1];
1564 const MFloat p1[nDim] = {wallDistSurfCoordsAll[0][globalWallId1], wallDistSurfCoordsAll[1][globalWallId1]};
1565 const MFloat p2[nDim] = {wallDistSurfCoordsAll[0][globalWallId2], wallDistSurfCoordsAll[1][globalWallId2]};
1566 const MFloat pTarget[nDim] = {m_cells->coordinates[0][cellId], m_cells->coordinates[1][cellId]};
1567 MFloat fraction;
1568 shortestDistanceToLineElement(p1, p2, pTarget, distance, fraction);
1569 // Sanity check
1570 if(fraction > 0.5 + 1e-8) {
1571 stringstream msg;
1572 msg << "ERROR: p1=" << p1[0] << "|" << p1[1] << " p2=" << p2[0] << "|" << p2[1] << " pTarget=" << pTarget[0]
1573 << "|" << pTarget[1] << " distance=" << distance << " fraction=" << fraction << endl;
1574 cout << msg.str();
1575 mTerm(1, "Fraction is >0.5, so nearest grid point is supposed to be p2 and not p1!");
1576 }
1577 weight = 1.0;
1578 } else { // grid point has neighbors to both sides (Case 1 & Case 2)
1579 // globalCellId1 lies in between the grid points globalWallId1 and globalWallId2
1580 // globalCellId2 lies in between the grid points globalWallId1 and globalWallId3 or
1581 // between the grid points globalWallId1_ and globalWallId3
1582 MInt globalWallId2;
1583 MInt globalWallId3;
1584 if(nfound == 1 /*&& isStartEndPointAll[globalWallId1]==0*/) {
1585 // Case 2:
1586 globalWallId2 = globalWallId1 - 1;
1587 globalWallId3 = globalWallId1 + 1;
1588 globalCellId1 = globalWallId2;
1589 globalCellId1 += cellIdShifts[globalWallId2];
1590 globalCellId2 = globalWallId1;
1591 globalCellId2 += cellIdShifts[globalWallId1];
1592 } else if(nfound == 2) { // globalWallId1 and globalWallId1_ coincide
1593 // Case 3:
1594 if(abs(isStartEndPointAll[globalWallId1]) + abs(isStartEndPointAll[globalWallId1_]) != 2) {
1595 // It can be a interface point which exists on two different solvers
1596 const MInt blockId1 = blockIds[getDomainId(globalWallId1)];
1597 const MInt blockId1_ = blockIds[getDomainId(globalWallId1_)];
1598 if(blockId1 == blockId1_) {
1599 cout << m_solver->domainId() << "|" << cellId << " " << globalWallId1 << " " << globalWallId1_ << " "
1600 << globalWallId1__ << ": cellCoord=" << setprecision(12) << m_cells->coordinates[0][cellId] << "|"
1601 << m_cells->coordinates[1][cellId] << "; wallCoord1=" << wallDistSurfCoordsAll[0][globalWallId1] << "|"
1602 << wallDistSurfCoordsAll[1][globalWallId1]
1603 << "; wallCoords1_=" << wallDistSurfCoordsAll[0][globalWallId1_] << "|"
1604 << wallDistSurfCoordsAll[1][globalWallId1_]
1605 << "; wallCoords1__=" << wallDistSurfCoordsAll[0][globalWallId1__] << "|"
1606 << wallDistSurfCoordsAll[1][globalWallId1__] << endl;
1607
1608 cout << "ERROR: (x|y)_1=" << setprecision(12) << wallDistSurfCoordsAll[0][globalWallId1] << "|"
1609 << wallDistSurfCoordsAll[1][globalWallId1] << " (x|y)_2=" << wallDistSurfCoordsAll[0][globalWallId1_]
1610 << "|" << wallDistSurfCoordsAll[1][globalWallId1_] << endl;
1611 mTerm(1, "One grid point found twice. Both grid points have neighbors to both sides, which is unexpected!");
1612 }
1613 if(blockId1 != myBlockId && blockId1_ != myBlockId) mTerm(1, "This case is not implemented yet!");
1614 if(blockId1_ == myBlockId) globalWallId1 = globalWallId1_;
1615 // Case 2:
1616 globalWallId2 = globalWallId1 - 1;
1617 globalWallId3 = globalWallId1 + 1;
1618 globalCellId1 = globalWallId2;
1619 globalCellId1 += cellIdShifts[globalWallId2];
1620 globalCellId2 = globalWallId1;
1621 globalCellId2 += cellIdShifts[globalWallId1];
1622 } else {
1623 globalWallId2 = globalWallId1 + isStartEndPointAll[globalWallId1];
1624 globalWallId3 = globalWallId1_ + isStartEndPointAll[globalWallId1_];
1625 globalCellId1 = (isStartEndPointAll[globalWallId1] == -1) ? globalWallId1 - 1 : globalWallId1;
1626 globalCellId1 += cellIdShifts[globalWallId1];
1627 globalCellId2 = (isStartEndPointAll[globalWallId1_] == -1) ? globalWallId1_ - 1 : globalWallId1_;
1628 globalCellId2 += cellIdShifts[globalWallId1_];
1629 }
1630 } else { // nfound==3
1631 // Here we could also check if the point belongs to different solvers
1632 stringstream msg;
1633 msg << setprecision(12) << "p1=" << wallDistSurfCoordsAll[0][globalWallId1] << "|"
1634 << wallDistSurfCoordsAll[1][globalWallId1] << ", p2=" << wallDistSurfCoordsAll[0][globalWallId1_] << "|"
1635 << wallDistSurfCoordsAll[1][globalWallId1_] << ", p3=" << wallDistSurfCoordsAll[0][globalWallId1__] << "|"
1636 << wallDistSurfCoordsAll[1][globalWallId1__] << endl;
1637 cout << msg.str();
1638 mTerm(1, "Three wall grid points coincide. This situation is unknown!");
1639 }
1640
1641 // Now compute distances to both line elements
1642 MFloat distances[2];
1643 MFloat fractions[2];
1644 MFloat lineLengths[2];
1645 const MFloat p1[nDim] = {wallDistSurfCoordsAll[0][globalWallId1], wallDistSurfCoordsAll[1][globalWallId1]};
1646 const MFloat p2[nDim] = {wallDistSurfCoordsAll[0][globalWallId2], wallDistSurfCoordsAll[1][globalWallId2]};
1647 const MFloat p3[nDim] = {wallDistSurfCoordsAll[0][globalWallId3], wallDistSurfCoordsAll[1][globalWallId3]};
1648 const MFloat pTarget[nDim] = {m_cells->coordinates[0][cellId], m_cells->coordinates[1][cellId]};
1649 lineLengths[0] = shortestDistanceToLineElement(p1, p2, pTarget, distances[0], fractions[0]);
1650 lineLengths[1] = shortestDistanceToLineElement(p1, p3, pTarget, distances[1], fractions[1]);
1651 const MInt i = distances[1] < distances[0];
1652 // Sanity check
1653 if(fractions[i] > 0.5 + 1e-8) {
1654 stringstream msg;
1655 msg << "Error p1=" << p1[0] << "|" << p1[1] << " p2=" << p2[0] << "|" << p2[1] << " p3=" << p3[0] << "|"
1656 << p3[1] << " pTarget=" << pTarget[0] << "|" << pTarget[1] << " distance=" << distance
1657 << " fraction=" << fractions[i] << " " << fractions[1 - i] << endl;
1658 cout << msg.str();
1659 mTerm(1, "Fraction >0.5, so nearest grid point can not be p1!");
1660 }
1661 // If wall is strongly concave, things might get nasty
1662 weight = (fractions[i] * lineLengths[i] + 0.5 * lineLengths[1 - i]) / (0.5 * (lineLengths[0] + lineLengths[1]));
1663 distance = distances[i];
1664 if(i == 1) {
1665 // Swap so that weight belongs to globalCellId1; weight of globalCellId2 is 1-weight
1666 const MInt temp = globalCellId1;
1667 globalCellId1 = globalCellId2;
1668 globalCellId2 = temp;
1669 }
1670 }
1671
1672 if(distance < r_cells_fq_distance[cellId]) {
1673 ASSERT(globalCellId1 != -1, "");
1674 // In case one globalId is -1, replace that id with the other one. Since weight of such a cell
1675 // is zero, it doesn't matter
1676 globalCellId2 = (globalCellId2 == -1) ? globalCellId1 : globalCellId2;
1677 r_cellId2globalWallId[cellId] = std::make_pair(globalCellId1, globalCellId2);
1678 r_weighting[cellId] = weight;
1679 r_cells_fq_distance[cellId] = distance;
1680
1681 // DEBUG
1682 // m_cells->fq[FQ->VAR1][cellId] = globalCellId1*weight + (1-weight)*globalCellId2;
1683 }
1684 }
1685
1686 // Convert displs based on grid points to displs based on cells
1687 for(MInt i = 1; i < (signed)r_displs.size(); ++i) {
1688 r_displs[i] = r_displs[i] + cellIdShifts[r_displs[i]];
1689 }
1690
1691 // Some debug output
1692#ifndef NDEBUG
1693 MInt noNearWallCells = 0;
1694 // globalWallIds: store the positions of all wall surfaces in wallDistSurfCoordsAll, which are relavant
1695 // for current rank
1696 std::set<MInt> globalWallIds;
1697 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1698 if(r_cellId2globalWallId[cellId].first != -1) {
1699 ++noNearWallCells;
1700 globalWallIds.insert(r_cellId2globalWallId[cellId].first);
1701 }
1702 if(r_cellId2globalWallId[cellId].second != -1) {
1703 globalWallIds.insert(r_cellId2globalWallId[cellId].second);
1704 }
1705 }
1706
1707 cout << "::computeDistance2Map: rank=" << m_solver->domainId() << " cellmemSize=" << cellmemSize
1708 << " noTotalWallPoints=" << noTotalWallPoints << " #cells near wall=" << noNearWallCells
1709 << " #globalWallIds=" << globalWallIds.size() << endl;
1710#endif
1711}
StructuredBndryCnd2D::shortestDistanceToLineElement
MFloat shortestDistanceToLineElement(const MFloat(&)[nDim], const MFloat(&)[nDim], const MFloat(&)[nDim], MFloat &, MFloat &)
Definition: fvstructuredbndrycnd2d.cpp:1100
id
MInt id
Definition: maiatypes.h:71
MPI_Allgatherv
int MPI_Allgatherv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, const int recvcounts[], const int displs[], MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgatherv
Definition: mpioverride.cpp:1315
MPI_Allgather
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
Definition: mpioverride.cpp:1274
maia::math::distance
MFloat distance(const MFloat *a, const MFloat *b)
Definition: maiamath.h:249
KDtree
Definition: kdtree.h:73
Point
Definition: pointbox.h:20
a
Definition: contexttypes.h:19

◆ computeFrictionCoef()

template<MBool isRans>
virtual void StructuredBndryCnd2D< isRans >::computeFrictionCoef ( )
inlineoverridevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 84 of file fvstructuredbndrycnd2d.h.

84{ computeFrictionPressureCoef_<false, true, false>(false, false); };

◆ computeFrictionPressureCoef()

template<MBool isRans>
virtual void StructuredBndryCnd2D< isRans >::computeFrictionPressureCoef ( MBool  computePower)
inlineoverridevirtual

Implements StructuredBndryCnd< 2 >.

Definition at line 81 of file fvstructuredbndrycnd2d.h.

81 {
82 computeFrictionPressureCoef_<true, true, true>(false, computePower);
83 }

◆ computeFrictionPressureCoef_()

template<MBool isRans>
template<MBool calcCp, MBool calcCf, MBool calcIntegrals>
template void StructuredBndryCnd2D< isRans >::computeFrictionPressureCoef_< true, true, true > ( const MBool  auxDataWindow = false,
const MBool  computePower = false 
)

Definition at line 2266 of file fvstructuredbndrycnd2d.cpp.

2266 {
2267 static_assert(calcCp || calcCf, "Something went wrong");
2268
2269 // References for convenience
2270 const MInt noForceCoefs = m_solver->m_noForceDataFields;
2271 auto& m_forceCoef = m_solver->m_forceCoef;
2272
2273 const MFloat UT = m_solver->m_Ma * sqrt(PV->TInfinity);
2274 const MFloat stagnationPressure = (CV->rhoInfinity * F1B2 * POW2(UT));
2275 const MFloat fstagnationEnergy = F1 / (CV->rhoInfinity * F1B2 * POW3(UT));
2276
2277 if(calcIntegrals) {
2278 // Only those for which auxData is requested in the property file; not those which are, e.g., required
2279 // because of k-epsilon rans model
2280 const MInt noWalls = m_solver->m_windowInfo->m_auxDataWindowIds.size();
2281
2282 // reset values to zero for the calculation
2283 for(MInt i = 0; i < (MInt)noForceCoefs * noWalls; ++i) {
2284 m_forceCoef[i] = F0;
2285 }
2286 }
2287
2288 // loop over all maps
2289 for(MInt map = 0; map < (MInt)m_auxDataMap.size(); ++map) {
2290 //
2291 if(auxDataWindows) {
2292 MBool takeIt = false;
2293 for(auto& m : m_solver->m_windowInfo->m_auxDataWindowIds) {
2294 if(m_auxDataMap[map]->Id2 == m.second) {
2295 takeIt = true;
2296 break;
2297 }
2298 }
2299 if(!takeIt) continue;
2300 }
2301
2302 const MBool interface = m_auxDataMap[map]->BC >= 6000 && m_auxDataMap[map]->BC < 6010;
2303 const auto calc_cp_cf_ = interface ? &StructuredBndryCnd2D<isRans>::calc_cp_cf<calcCp, calcCf, true>
2304 : &StructuredBndryCnd2D<isRans>::calc_cp_cf<calcCp, calcCf, false>;
2305
2306 const MInt mapOffsetCf = m_cells->cfOffsets[map];
2307 const MInt mapOffsetCp = m_cells->cpOffsets[map];
2308 MInt* start = m_auxDataMap[map]->start1;
2309 MInt* end = m_auxDataMap[map]->end1;
2310
2311 MInt n = 0;
2312 MInt normalDir = 0;
2313 // area
2314 MFloat area = F0;
2315 // skin-friction
2316 MFloat cf[nDim] = {};
2317 // pressure coefficient
2318 MFloat cp[nDim] = {};
2319
2320 MFloat powerp[2] = {F0, F0};
2321 MFloat powerv[2] = {F0, F0};
2322
2323 // indices
2324 MInt n10[nDim] = {};
2325 MInt n01[nDim] = {};
2326 MInt n1m1[nDim] = {};
2327 MInt pCoordDir[(2 * (nDim - 1) - 1) * nDim] = {};
2328 MInt i = 0, j = 0;
2329
2330 // output for cp & cf-vector
2331 MFloat output[calcCp + nDim * calcCf];
2332
2333 // determine the wall
2334 if((m_auxDataMap[map]->face) % 2 == 0) {
2335 n = -1; // bottom wall
2336 normalDir = m_auxDataMap[map]->face / 2;
2337 i = start[normalDir];
2338 } else {
2339 n = 1; // top wall
2340 normalDir = (m_auxDataMap[map]->face - 1) / 2; // TODO_SS labels:FV Check if this and
2341 i = end[normalDir] - 1; // the -1 at this line is necessary
2342 }
2343 // determine the normals
2344 n1m1[normalDir] = -1 * n;
2345 n10[normalDir] = (MInt)(0.5 + (0.5 * (MFloat)n));
2346 n01[normalDir] = (MInt)(-0.5 + (0.5 * (MFloat)n));
2347
2348 const MInt firstTangential = 1 - normalDir;
2349 // index shift for other surface point
2350 pCoordDir[firstTangential] = 1;
2351 // size in tangential direction
2352 const MInt sizeJ = end[firstTangential] - start[firstTangential];
2353 MInt& reali = (normalDir == 0) ? i : j;
2354 MInt& realj = (normalDir == 0) ? j : i;
2355
2356 // loop over the tangential direction of the wall
2357 MInt jj = 0;
2358 for(j = start[firstTangential]; j < end[firstTangential]; j++) {
2359 // get id of boundary cell and cell above that one for extrapolation
2360 const MInt cellId = cellIndex(reali, realj);
2361 const MInt cellIdP1 = cellIndex(reali + n1m1[0], realj + n1m1[1]);
2362
2363 // get point id and ids of three other surface corners
2364 const MInt pIJ = getPointIdFromCell(reali + n10[0], realj + n10[1]);
2365 const MInt pIJP = getPointIdFromPoint(pIJ, pCoordDir[0], pCoordDir[1]);
2366
2367 // Actual calculation of cp & cf
2368 (this->*calc_cp_cf_)(cellId, cellIdP1, pIJ, pIJP, output);
2369
2370 if(calcCp)
2371 // save pressure to map
2372 m_cells->cp[mapOffsetCp + jj] = output[0]; // cpn;
2373
2374 if(calcCf) {
2375 // save skin-friction to map
2376 m_cells->cf[mapOffsetCf + 0 * sizeJ + jj] = output[calcCp + 0]; // cfx;
2377 m_cells->cf[mapOffsetCf + 1 * sizeJ + jj] = output[calcCp + 1]; // cfy;
2378 }
2379
2380 if((calcCf && m_solver->m_detailAuxData) || calcIntegrals) {
2381 // first get the position of the wall (reference!!)
2382 // x is always used as coordinate in normal direction
2383 const MFloat xRef = 0.5 * (m_grid->m_coordinates[0][pIJ] + m_grid->m_coordinates[0][pIJP]);
2384 const MFloat yRef = 0.5 * (m_grid->m_coordinates[1][pIJ] + m_grid->m_coordinates[1][pIJP]);
2385
2386 // this should be the valid method as in TFS
2387 // cp:
2388 // sum up all lengths and all cps with their surface width contribution
2389 const MFloat dxidx = m_cells->surfaceMetrics[nDim * normalDir + 0][cellIndex(reali + n01[0], realj + n01[1])];
2390 const MFloat dxidy = m_cells->surfaceMetrics[nDim * normalDir + 1][cellIndex(reali + n01[0], realj + n01[1])];
2391
2392 if(calcCf && m_solver->m_detailAuxData) {
2393 // save surface contributions
2394 m_cells->cf[mapOffsetCf + 2 * sizeJ + jj] = dxidx;
2395 m_cells->cf[mapOffsetCf + 3 * sizeJ + jj] = dxidy;
2396
2397 // save surface coordinates
2398 m_cells->cf[mapOffsetCf + 4 * sizeJ + jj] = xRef;
2399 m_cells->cf[mapOffsetCf + 5 * sizeJ + jj] = yRef;
2400 }
2401
2402 if(calcIntegrals) {
2403 MFloat considerValue = F1;
2404 if(m_solver->m_auxDataCoordinateLimits) {
2405 if(m_solver->m_auxDataLimits[0] <= xRef && xRef <= m_solver->m_auxDataLimits[1]
2406 && m_solver->m_auxDataLimits[2] <= yRef && yRef <= m_solver->m_auxDataLimits[3]) {
2407 considerValue = F1;
2408 } else {
2409 considerValue = F0;
2410 }
2411 }
2412
2413 if(calcCp) {
2414 cp[0] += (-1.0) * output[0] * dxidx * considerValue;
2415 cp[1] += (-1.0) * output[0] * dxidy * considerValue;
2416 }
2417 if(calcCf) {
2418 // cf
2419 cf[0] += output[calcCp + 0] * (sqrt(dxidx * dxidx + dxidy * dxidy)) * considerValue;
2420 cf[1] += output[calcCp + 1] * (sqrt(dxidx * dxidx + dxidy * dxidy)) * considerValue;
2421 }
2422
2423 if(computePower) {
2424 const MFloat uWall =
2425 m_solver->m_movingGrid ? F1B2 * (m_grid->m_velocity[0][pIJ] + m_grid->m_velocity[0][pIJP]) : F0;
2426 const MFloat vWall =
2427 m_solver->m_movingGrid ? F1B2 * (m_grid->m_velocity[1][pIJ] + m_grid->m_velocity[1][pIJP]) : F0;
2428
2429 const MFloat dp = output[0] * stagnationPressure;
2430
2431 const MFloat P_px = (-1.0) * dp * uWall * fstagnationEnergy;
2432 const MFloat P_py = (-1.0) * dp * vWall * fstagnationEnergy;
2433
2434 const MFloat taux = output[calcCp + 0] * stagnationPressure;
2435 const MFloat tauy = output[calcCp + 1] * stagnationPressure;
2436
2437 const MFloat P_cfx = (taux)*uWall * fstagnationEnergy;
2438 const MFloat P_cfy = (tauy)*vWall * fstagnationEnergy;
2439
2440 powerp[0] += P_px * dxidx * considerValue;
2441 powerp[1] += P_py * dxidy * considerValue;
2442
2443 powerv[0] += P_cfx * (sqrt(dxidx * dxidx + dxidy * dxidy)) * considerValue;
2444 powerv[1] += P_cfy * (sqrt(dxidx * dxidx + dxidy * dxidy)) * considerValue;
2445 }
2446 // area
2447 area += sqrt(dxidx * dxidx + dxidy * dxidy) * considerValue;
2448 }
2449 }
2450 ++jj;
2451 }
2452
2453 if(calcIntegrals) {
2454 // now add to lift and drag coefficients
2455 MInt count = 0;
2456 for(auto it = m_solver->m_windowInfo->m_auxDataWindowIds.cbegin();
2457 it != m_solver->m_windowInfo->m_auxDataWindowIds.cend();
2458 ++it) {
2459 if(m_auxDataMap[map]->Id2 == it->second) {
2460 m_forceCoef[count * noForceCoefs + 0] += cf[0];
2461 m_forceCoef[count * noForceCoefs + 1] += cf[1];
2462 m_forceCoef[count * noForceCoefs + 2] += 0.0; // the compressible part of the stress tensor - x
2463 m_forceCoef[count * noForceCoefs + 3] += 0.0; // the compressible part of the stress tensor - y
2464 m_forceCoef[count * noForceCoefs + 4] += cp[0];
2465 m_forceCoef[count * noForceCoefs + 5] += cp[1];
2466 m_forceCoef[count * noForceCoefs + 6] += area;
2467 }
2468
2469 if(computePower) {
2470 m_forceCoef[count * noForceCoefs + 7] += powerv[0];
2471 m_forceCoef[count * noForceCoefs + 8] += powerv[1];
2472 m_forceCoef[count * noForceCoefs + 9] += powerp[0];
2473 m_forceCoef[count * noForceCoefs + 10] += powerp[1];
2474 }
2475
2476 count++;
2477 }
2478 }
2479 }
2480}
StructuredBndryCnd2D
Class for the 2D stuctured boundary conditions.
Definition: fvstructuredbndrycnd2d.h:21
StructuredBndryCnd< 2 >::m_auxDataMap
std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > & m_auxDataMap
Definition: fvstructuredbndrycnd.h:63
POW3
constexpr Real POW3(const Real x)
Definition: functions.h:123

◆ computeLocalExtendedDistancesAndSetComm()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::computeLocalExtendedDistancesAndSetComm
overridevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 1157 of file fvstructuredbndrycnd2d.cpp.

1157 {
1158 TRACE();
1159
1160 // Store solid maps and fluid-porous maps in seperate lists
1161 vector<unique_ptr<StructuredWindowMap<nDim>>> wallDistInfo;
1162 vector<unique_ptr<StructuredWindowMap<nDim>>> FPDistInfo;
1163 for(auto& auxDataMap : m_auxDataMap) {
1164 if(auxDataMap->isFluidPorousInterface) {
1165 unique_ptr<StructuredWindowMap<nDim>> temp = make_unique<StructuredWindowMap<nDim>>();
1166 m_solver->m_windowInfo->mapCpy(auxDataMap, temp);
1167 FPDistInfo.push_back(move(temp));
1168 } else if((int)(auxDataMap->BC / 1000.0) == 1) {
1169 unique_ptr<StructuredWindowMap<nDim>> temp = make_unique<StructuredWindowMap<nDim>>();
1170 m_solver->m_windowInfo->mapCpy(auxDataMap, temp);
1171 wallDistInfo.push_back(move(auxDataMap));
1172 }
1173 }
1174#ifndef NDEBUG
1175 std::cout << "DomainId=" << m_solver->domainId() << ": " << FPDistInfo.size() << " FP maps & " << wallDistInfo.size()
1176 << " wall maps." << endl;
1177#endif
1178
1179 //
1180 MInt noRanks;
1181 MPI_Comm_size(m_solver->m_StructuredComm, &noRanks);
1182
1183 // Allocate temporary space; the distance to solid wall is temporarily saved in m_cells->fq
1184 ScratchSpace<MFloat> FP_distance(m_noCells, AT_, "FP_distance");
1185
1186 // Solid
1187 m_wallSendcounts.assign(noRanks, 0);
1188 // r_cellId2globalWallId: mapping of cellId to the positions of the two cells which encloses the nearest point at the
1189 // wall
1190 std::vector<std::pair<MInt, MInt>> cellId2globalWallId; //(m_noCells, std::make_pair(-1, -1));
1191 std::vector<MInt> wallDispls;
1192 // weighting: weighting of the wall cell which is closest to a given cell; the weighting of its neighbor is
1193 // 1-weighting
1194 std::vector<MFloat> wallWeighting; //(m_noCells, 0.0);
1195 computeDistance2Map(wallDistInfo, m_cells->fq[FQ->WALLDISTANCE], cellId2globalWallId, wallDispls, wallWeighting);
1196
1197 // Porous
1198 m_FPSendcounts.assign(noRanks, 0);
1199 // r_cellId2globalWallId: mapping of cellId to the positions of the two cells which encloses the nearest point at the
1200 // wall
1201 std::vector<std::pair<MInt, MInt>> cellId2globalFPId; //(m_noCells, std::make_pair(-1,-1));
1202 std::vector<MInt> FPDispls;
1203 // weighting: weighting of the wall cell which is closest to a given cell; the weighting of its neighbor is
1204 // 1-weighting
1205 std::vector<MFloat> FPWeighting; //(m_noCells, 0.0);
1206 computeDistance2Map(FPDistInfo, &FP_distance[0], cellId2globalFPId, FPDispls, FPWeighting);
1207
1209 // Before we reset we reset the communication information of porous solvers about fluid-porous interfaces,
1210 // we save that information
1211 std::vector<std::pair<MInt, MInt>> cellId2globalFPId_(cellId2globalFPId.begin(), cellId2globalFPId.end());
1212 if(m_solver->m_blockType != "porous") {
1213 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1214 cellId2globalFPId_[cellId] = std::make_pair(-1, -1);
1215 }
1216 }
1217 setUpNearMapComm(cellId2globalFPId_, FPDispls, FPWeighting, m_FPCellId2recvCell_porous, m_FPSendcounts_porous,
1218 m_FPSendCells_porous);
1220
1221 // Here we need to modify the walldistance; if we are inside porous solver set wall FPdistance to max value
1222 modifyFPDistance(FPDistInfo, &FP_distance[0], cellId2globalFPId, wallWeighting);
1223
1224 //
1225 getCloserMap(m_cells->fq[FQ->WALLDISTANCE], cellId2globalWallId, &FP_distance[0], cellId2globalFPId,
1226 m_cells->fq[FQ->WALLDISTANCE]);
1227
1228 // Solid
1229 setUpNearMapComm(cellId2globalWallId, wallDispls, wallWeighting, m_wallCellId2recvCell, m_wallSendcounts,
1230 m_wallSendCells);
1231
1232 // Porous
1233 setUpNearMapComm(cellId2globalFPId, FPDispls, FPWeighting, m_FPCellId2recvCell, m_FPSendcounts, m_FPSendCells);
1234}
StructuredBndryCnd2D::getCloserMap
void getCloserMap(const MFloat *const, std::vector< std::pair< MInt, MInt > > &, const MFloat *const, std::vector< std::pair< MInt, MInt > > &, MFloat *const, T comparator={})
Definition: fvstructuredbndrycnd2d.cpp:1876
StructuredBndryCnd2D::setUpNearMapComm
void setUpNearMapComm(const std::vector< std::pair< MInt, MInt > > &, const std::vector< MInt > &, const std::vector< MFloat > &, std::map< MInt, std::tuple< MInt, MInt, MFloat > > &, std::vector< MInt > &, std::vector< MInt > &)
Definition: fvstructuredbndrycnd2d.cpp:1914
StructuredBndryCnd2D::computeDistance2Map
void computeDistance2Map(const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, MFloat *const, std::vector< std::pair< MInt, MInt > > &, std::vector< MInt > &, std::vector< MFloat > &)
Compute shortest distance to given set of maps.
Definition: fvstructuredbndrycnd2d.cpp:1366
StructuredBndryCnd2D::modifyFPDistance
void modifyFPDistance(const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, MFloat *const, std::vector< std::pair< MInt, MInt > > &, const std::vector< MFloat > &)
Modifies the fluid-porous distance.
Definition: fvstructuredbndrycnd2d.cpp:1244
StructuredBndryCnd< 2 >::m_FPSendcounts_porous
std::vector< MInt > m_FPSendcounts_porous
Definition: fvstructuredbndrycnd.h:206
StructuredBndryCnd< 2 >::m_wallCellId2recvCell
std::map< MInt, std::tuple< MInt, MInt, MFloat > > m_wallCellId2recvCell
Definition: fvstructuredbndrycnd.h:199
StructuredBndryCnd< 2 >::m_wallSendcounts
std::vector< MInt > m_wallSendcounts
Definition: fvstructuredbndrycnd.h:198
StructuredBndryCnd< 2 >::m_FPCellId2recvCell_porous
std::map< MInt, std::tuple< MInt, MInt, MFloat > > m_FPCellId2recvCell_porous
Definition: fvstructuredbndrycnd.h:207
StructuredBndryCnd< 2 >::m_FPSendCells
std::vector< MInt > m_FPSendCells
Definition: fvstructuredbndrycnd.h:201
StructuredBndryCnd< 2 >::m_wallSendCells
std::vector< MInt > m_wallSendCells
Definition: fvstructuredbndrycnd.h:197
StructuredBndryCnd< 2 >::m_FPSendcounts
std::vector< MInt > m_FPSendcounts
Definition: fvstructuredbndrycnd.h:202
StructuredBndryCnd< 2 >::m_FPCellId2recvCell
std::map< MInt, std::tuple< MInt, MInt, MFloat > > m_FPCellId2recvCell
Definition: fvstructuredbndrycnd.h:203
StructuredBndryCnd< 2 >::m_FPSendCells_porous
std::vector< MInt > m_FPSendCells_porous
Definition: fvstructuredbndrycnd.h:205

◆ computeLocalWallDistances()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::computeLocalWallDistances
overridevirtual

Set distance, nearest neigbor cell identifiers & weighting

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 741 of file fvstructuredbndrycnd2d.cpp.

741 {
742 // Initialize array with high numbers
743 MFloat maxWallDistance = numeric_limits<MFloat>::max();
744 maxWallDistance = Context::getSolverProperty<MFloat>("maxWallDistance", m_solverId, AT_, &maxWallDistance);
745 for(MInt id = 0; id < m_noCells; ++id) {
746 m_cells->fq[FQ->WALLDISTANCE][id] = maxWallDistance;
747 }
748
749 const vector<unique_ptr<StructuredWindowMap<nDim>>>& wallDistInfo =
750 m_auxDataMap; // m_solver->m_windowInfo->m_wallDistInfoMap;
751 const MInt noWallDistInfo = wallDistInfo.size(); // contains the size of the maps
752 std::vector<MInt> wallNormalIndex(noWallDistInfo);
753 std::vector<MInt> normalDir(noWallDistInfo);
754 std::vector<std::array<MInt, 2 * nDim>> startEndTangential(noWallDistInfo);
755 // Number of grid points on current rank to store
756 MInt cellmemSize = 0;
757 for(MInt nn = 0; nn < noWallDistInfo; ++nn) {
758 // start & end represent exactly the range of grid point indices, which form the map/boundary
759 const MInt* const start = wallDistInfo[nn]->start1;
760 const MInt* const end = wallDistInfo[nn]->end1;
761 const MInt face = wallDistInfo[nn]->face;
762 normalDir[nn] = face / 2;
763 wallNormalIndex[nn] = start[normalDir[nn]];
764
765 MInt size = 1;
766 for(MInt dim = 0; dim < nDim - 1; ++dim) {
767 const MInt tangentialDir = (normalDir[nn] + (dim + 1)) % nDim;
768 startEndTangential[nn][dim * 2 + 0] = start[tangentialDir]; // + m_noGhostLayers;
769 startEndTangential[nn][dim * 2 + 1] =
770 end[tangentialDir] + 1; //+1, because grid coords not cell coords// - m_noGhostLayers;
771 size *= startEndTangential[nn][dim * 2 + 1] - startEndTangential[nn][dim * 2 + 0];
772 }
773 cellmemSize += size;
774#ifndef NDEBUG
775 cout << "Debug output: wallDistInfo=" << nn << ": face=" << face << "; normalDir=" << normalDir[nn]
776 << "; wallNormalIndex=" << wallNormalIndex[nn] << "; startEndTangential=" << startEndTangential[nn][0] << "|"
777 << startEndTangential[nn][1] << "normal start/end=" << start[normalDir[nn]] << "|" << end[normalDir[nn]]
778 << endl;
779#endif
780 }
781
782 // 2.2) allocate the space for all the wall grid coordinates in Scratch!
783 // memory will not be needed later!
784 MFloatScratchSpace coordGridMem(std::max(1, nDim * cellmemSize), AT_, "coordGridMem");
785 MFloatPointerScratchSpace wallDistSurfCoords(nDim, AT_, "wallDistSurfCoords");
786 for(MInt dim = 0; dim < nDim; ++dim)
787 wallDistSurfCoords[dim] = &coordGridMem[dim * cellmemSize];
788 // Store if a grid point has a neighbor to the left or right on the same domain
789 MIntScratchSpace isStartEndPoint(std::max(1, cellmemSize), AT_, "isStartEndPoint");
790
791 // 4) Put wall grid points and information about existence of neighbor grid point in temporary buffers;
792 MInt cnt = 0;
793 for(MInt nn = 0; nn < noWallDistInfo; ++nn) {
794 MInt tIdx;
795 MInt& i = (normalDir[nn] == 0) ? wallNormalIndex[nn] : tIdx;
796 MInt& j = (normalDir[nn] == 0) ? tIdx : wallNormalIndex[nn];
797 for(tIdx = startEndTangential[nn][0]; tIdx < startEndTangential[nn][1]; ++tIdx) {
798 const MInt p = getPointIdFromCell(i, j);
799 wallDistSurfCoords[0][cnt] = m_grid->m_coordinates[0][p];
800 wallDistSurfCoords[1][cnt] = m_grid->m_coordinates[1][p];
801 // The case no neighbor to the left and right can not occur
802 if(tIdx == startEndTangential[nn][0])
803 isStartEndPoint[cnt] = 1; // no neighbor to the left on this rank
804 else if(tIdx == startEndTangential[nn][1] - 1)
805 isStartEndPoint[cnt] = -1; // no neighbor to the right on this rank
806 else
807 isStartEndPoint[cnt] =
808 0; // TODO_SS labels:FV,GRID,totest is this necessary, or will it be default initialized anyway?
809 ++cnt;
810 }
811 }
812
813 // Broadcast all wall grid coordinates and isStartEndPoint to all ranks
814 MInt noRanks;
815 MPI_Comm_size(m_solver->m_StructuredComm, &noRanks);
816 // recvcounts: How many wall grid coordinates to receive from each domain
817 std::vector<MInt> recvcounts(noRanks);
818 MPI_Allgather(&cellmemSize, 1, MPI_INT, recvcounts.data(), 1, MPI_INT, m_solver->m_StructuredComm, AT_, "cellmemSize",
819 "recvcounts");
820 std::vector<MInt> displs(noRanks);
821 for(MInt r = 0; r < noRanks - 1; ++r) {
822 displs[r + 1] = displs[r] + recvcounts[r];
823 }
824 const MInt noTotalWallPoints = displs[noRanks - 1] + recvcounts[noRanks - 1];
825 // wallDistSurfCoordsAll: contains the wall grid coordinates of all ranks
826 std::vector<std::vector<MFloat>> wallDistSurfCoordsAll(nDim);
827 for(MInt dim = 0; dim < nDim; ++dim) {
828 wallDistSurfCoordsAll[dim].resize(noTotalWallPoints);
829 MPI_Allgatherv(&wallDistSurfCoords[dim][0], cellmemSize, MPI_DOUBLE, wallDistSurfCoordsAll[dim].data(),
830 recvcounts.data(), displs.data(), MPI_DOUBLE, m_solver->m_StructuredComm, AT_, "wallDistSurfCoords",
831 "wallDistSurfCoordsAll");
832 }
833 std::vector<MInt> isStartEndPointAll(noTotalWallPoints);
834 MPI_Allgatherv(&isStartEndPoint[0], cellmemSize, MPI_INT, isStartEndPointAll.data(), recvcounts.data(), displs.data(),
835 MPI_INT, m_solver->m_StructuredComm, AT_, "isStartEndPoint", "isStartEndPointAll");
836
837 // --> Now all ranks have all wall grid coordinates and the information if that grid point has any
838 // neighbor to the left or right on same domain
839
840 // globalWallIds: store the positions of all wall grid points in wallDistSurfCoordsAll, which are relevant
841 // for current rank
842 std::set<MInt> globalWallIds;
843 // cellId2globalWallId: mapping of cellId to the positions of the two cells which encloses the nearest point at the
844 // wall
845 std::vector<std::pair<MInt, MInt>> cellId2globalWallId(m_noCells, std::make_pair(-1, -1));
846 // weighting: weighting of the wall cell which is closest to a given cell; the weighting of its neighbor is
847 // 1-weighting
848 std::vector<MFloat> weighting(m_noCells, 0.0);
849
850 // Build k-d-tree for a quick search
851 // In the following Point<3> is used, but 3rd dimension is just dummy
852 std::vector<Point<3>> pts;
853 // cellIdShifts is used to retrieve the correct cellIds from the grid ids
854 std::vector<MInt> cellIdShifts(noTotalWallPoints + 1, 0);
855 for(MInt globalWallId = 0; globalWallId < noTotalWallPoints; ++globalWallId) {
856 Point<3> a(wallDistSurfCoordsAll[0][globalWallId], wallDistSurfCoordsAll[1][globalWallId], 0, globalWallId);
857 pts.push_back(a);
858 cellIdShifts[globalWallId + 1] = cellIdShifts[globalWallId] - 1 * (isStartEndPointAll[globalWallId] == -1);
859 }
860 // build up the tree
861 KDtree<3> tree(pts);
862 static constexpr MFloat radius = std::numeric_limits<MFloat>::max();
863 static constexpr MFloat eps = 1e-16; // I hope this is a suitable value
864
865 // go through all the cells an determine the closest distance
866 MInt results[3];
867 MFloat dists[3];
868 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
869 Point<3> pt(m_cells->coordinates[0][cellId], m_cells->coordinates[1][cellId], 0);
870 MBool overflow = false;
871 MInt nfound = tree.locatenearest(pt, radius, results, dists, 3, overflow);
872 if(nfound != 3) mTerm(1, "Come on, your mesh should have at least 3 wall grid points!");
873 const MInt globalWallId1 = results[0];
874 MInt globalWallId1_ = results[1];
875 MInt globalWallId1__ = results[2];
876 nfound = 1;
877 const MFloat r2 = POW2(wallDistSurfCoordsAll[0][globalWallId1] - wallDistSurfCoordsAll[0][globalWallId1_])
878 + POW2(wallDistSurfCoordsAll[1][globalWallId1] - wallDistSurfCoordsAll[1][globalWallId1_]);
879 const MFloat r3 = POW2(wallDistSurfCoordsAll[0][globalWallId1] - wallDistSurfCoordsAll[0][globalWallId1__])
880 + POW2(wallDistSurfCoordsAll[1][globalWallId1] - wallDistSurfCoordsAll[1][globalWallId1__]);
881 nfound += r2 < eps;
882 nfound += r3 < eps;
883 // It might happen that all three points have same distance to target point, but point 1 and 3 coincide and not
884 // 1 and 2
885 if(nfound == 2 && r3 < r2) {
886 const MInt temp = globalWallId1_;
887 globalWallId1_ = globalWallId1__;
888 globalWallId1__ = temp;
889 }
890
894 // Case 1: wall grid point has only one neighbor, because it is the begin or end of a solid wall
895 // Case 2: wall grid point has two neighbors, both on one rank (nfound==1)
896 // Case 3: wall grid point has two neighbors, but one neighbor is on a different rank (nfound==2)
897 MInt globalCellId1 = -1;
898 MInt globalCellId2 = -1;
899 MFloat distance{}, weight;
900 if(nfound == 1 && abs(isStartEndPointAll[globalWallId1]) == 1) {
901 // Case 1:
902 const MInt globalWallId2 = globalWallId1 + isStartEndPointAll[globalWallId1];
903 globalCellId1 = (isStartEndPointAll[globalWallId1] == -1) ? globalWallId1 - 1 : globalWallId1;
904 globalCellId1 += cellIdShifts[globalWallId1];
905 const MFloat p1[nDim] = {wallDistSurfCoordsAll[0][globalWallId1], wallDistSurfCoordsAll[1][globalWallId1]};
906 const MFloat p2[nDim] = {wallDistSurfCoordsAll[0][globalWallId2], wallDistSurfCoordsAll[1][globalWallId2]};
907 const MFloat pTarget[nDim] = {m_cells->coordinates[0][cellId], m_cells->coordinates[1][cellId]};
908 MFloat fraction;
909 shortestDistanceToLineElement(p1, p2, pTarget, distance, fraction);
910 // Sanity check
911 if(fraction > 0.5 + 1e-8) {
912 stringstream msg;
913 msg << "ERROR: p1=" << p1[0] << "|" << p1[1] << " p2=" << p2[0] << "|" << p2[1] << " pTarget=" << pTarget[0]
914 << "|" << pTarget[1] << " distance=" << distance << " fraction=" << fraction << endl;
915 cout << msg.str();
916 mTerm(1, "Fraction is >0.5, so nearest grid point is supposed to be p2 and not p1!");
917 }
918 weight = 1.0;
919 } else { // grid point has neighbors to both sides (Case 1 & Case 2)
920 // globalCellId1 lies in between the grid points globalWallId1 and globalWallId2
921 // globalCellId2 lies in between the grid points globalWallId1 and globalWallId3 or
922 // between the grid points globalWallId1_ and globalWallId3
923 MInt globalWallId2;
924 MInt globalWallId3;
925 if(nfound == 1) {
926 // Case 2:
927 globalWallId2 = globalWallId1 - 1;
928 globalWallId3 = globalWallId1 + 1;
929 globalCellId1 = globalWallId2;
930 globalCellId1 += cellIdShifts[globalWallId2];
931 globalCellId2 = globalWallId1;
932 globalCellId2 += cellIdShifts[globalWallId1];
933 } else if(nfound == 2) {
934 // Case 3:
935 if(abs(isStartEndPointAll[globalWallId1]) + abs(isStartEndPointAll[globalWallId1_]) != 2) {
936 cout << "ERROR: (x|y)_1=" << setprecision(12) << wallDistSurfCoordsAll[0][globalWallId1] << "|"
937 << wallDistSurfCoordsAll[1][globalWallId1] << " (x|y)_2=" << wallDistSurfCoordsAll[0][globalWallId1_]
938 << "|" << wallDistSurfCoordsAll[1][globalWallId1_] << endl;
939 mTerm(1, "One grid point found twice. Both grid points have neighbors to both sides, which is unexpected!");
940 }
941 globalWallId2 = globalWallId1 + isStartEndPointAll[globalWallId1];
942 globalWallId3 = globalWallId1_ + isStartEndPointAll[globalWallId1_];
943 globalCellId1 = (isStartEndPointAll[globalWallId1] == -1) ? globalWallId1 - 1 : globalWallId1;
944 globalCellId1 += cellIdShifts[globalWallId1];
945 globalCellId2 = (isStartEndPointAll[globalWallId1_] == -1) ? globalWallId1_ - 1 : globalWallId1_;
946 globalCellId2 += cellIdShifts[globalWallId1_];
947 } else { // nfound==3
948 stringstream msg;
949 msg << "p1=" << wallDistSurfCoordsAll[0][globalWallId1] << "|" << wallDistSurfCoordsAll[1][globalWallId1]
950 << ", p2=" << wallDistSurfCoordsAll[0][globalWallId1_] << "|" << wallDistSurfCoordsAll[1][globalWallId1_]
951 << ", p3=" << wallDistSurfCoordsAll[0][globalWallId1__] << "|" << wallDistSurfCoordsAll[1][globalWallId1__]
952 << endl;
953 cout << msg.str();
954 mTerm(1, "Three wall grid points coincide. This situation is unknown!");
955 }
956
957 // Now compute distances to both line elements
958 MFloat distances[2];
959 MFloat fractions[2];
960 MFloat lineLengths[2];
961 const MFloat p1[nDim] = {wallDistSurfCoordsAll[0][globalWallId1], wallDistSurfCoordsAll[1][globalWallId1]};
962 const MFloat p2[nDim] = {wallDistSurfCoordsAll[0][globalWallId2], wallDistSurfCoordsAll[1][globalWallId2]};
963 const MFloat p3[nDim] = {wallDistSurfCoordsAll[0][globalWallId3], wallDistSurfCoordsAll[1][globalWallId3]};
964 const MFloat pTarget[nDim] = {m_cells->coordinates[0][cellId], m_cells->coordinates[1][cellId]};
965 lineLengths[0] = shortestDistanceToLineElement(p1, p2, pTarget, distances[0], fractions[0]);
966 lineLengths[1] = shortestDistanceToLineElement(p1, p3, pTarget, distances[1], fractions[1]);
967 const MInt i = distances[1] < distances[0];
968 // Sanity check
969 if(fractions[i] > 0.5 + 1e-8) {
970 stringstream msg;
971 msg << "Error p1=" << p1[0] << "|" << p1[1] << " p2=" << p2[0] << "|" << p2[1] << " p3=" << p3[0] << "|"
972 << p3[1] << " pTarget=" << pTarget[0] << "|" << pTarget[1] << " distance=" << distance
973 << " fraction=" << fractions[i] << " " << fractions[1 - i] << endl;
974 cout << msg.str();
975 mTerm(1, "Fraction >0.5, so nearest grid point can not be p1!");
976 }
977 // If wall is strongly concave, things might get nasty
978 weight = (fractions[i] * lineLengths[i] + 0.5 * lineLengths[1 - i]) / (0.5 * (lineLengths[0] + lineLengths[1]));
979 distance = distances[i];
980 if(i == 1) {
981 // Swap so that weight belongs to globalCellId1; weight of globalCellId2 is 1-weight
982 const MInt temp = globalCellId1;
983 globalCellId1 = globalCellId2;
984 globalCellId2 = temp;
985 }
986 }
987
988 if(distance < m_cells->fq[FQ->WALLDISTANCE][cellId]) {
989 globalWallIds.insert(globalCellId1);
990 globalWallIds.insert(globalCellId2);
991 cellId2globalWallId[cellId] = std::make_pair(globalCellId1, globalCellId2);
992 weighting[cellId] = weight;
993 m_cells->fq[FQ->WALLDISTANCE][cellId] = distance;
994
995 // DEBUG
996 // m_cells->fq[FQ->VAR1][cellId] = globalCellId1*weight + (1-weight)*globalCellId2;
997 }
998 }
999 globalWallIds.erase(-1);
1000
1001 // Convert displs based on grid points to displs based on cells
1002 for(MInt i = 1; i < (signed)displs.size(); ++i) {
1003 displs[i] = displs[i] + cellIdShifts[displs[i]];
1004 }
1005 auto getDomainId = [&displs, noRanks](const MInt id, const MInt hint = 0) {
1006 for(MInt rank = hint; rank < (signed)displs.size() - 1; ++rank) {
1007 if(displs[rank + 1] > id) return rank;
1008 }
1009 return noRanks - 1;
1010 };
1011
1012 // globalWallId2localId: mapping of globalWallId to local wallId on current rank
1013 std::vector<MInt> globalWallId2localId((globalWallIds.empty()) ? 0 : *globalWallIds.rbegin() + 1);
1014 // recvWallIds: local wallIds from each rank, requested by current rank
1015 std::vector<MInt> recvWallIds(globalWallIds.size());
1016 std::vector<MInt> sendcounts(noRanks);
1017 MInt domainId_temp = 0;
1018 MInt localWallId = 0;
1019 for(auto it = globalWallIds.begin(); it != globalWallIds.end(); ++it, ++localWallId) {
1020 globalWallId2localId[*it] = localWallId;
1021 domainId_temp = getDomainId(*it, domainId_temp);
1022 // transform globalWallId "*it" to local wallId of rank domainId_temp
1023 recvWallIds[localWallId] = *it - displs[domainId_temp];
1024 ++sendcounts[domainId_temp];
1025 }
1026
1027 // m_wallCellId2recvCell: mapping from local cellId to the position, where its closest wall surface information is
1028 // stored
1029 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1030 if(cellId2globalWallId[cellId].first == -1 && cellId2globalWallId[cellId].second == -1) continue;
1031 // In case one globalWallId is -1, replace that id with the other one. Since weight of such a cell
1032 // is zero, it doesn't matter
1033 const MInt temp = (cellId2globalWallId[cellId].second == -1) ? cellId2globalWallId[cellId].first
1034 : cellId2globalWallId[cellId].second;
1035 auto temp2 = std::make_tuple(globalWallId2localId[cellId2globalWallId[cellId].first], globalWallId2localId[temp],
1036 weighting[cellId]);
1037 m_wallCellId2recvCell.insert({cellId, temp2});
1038 }
1039
1040 // Broadcast send/receive infos
1041 std::vector<MInt> snghbrs(noRanks);
1042 std::iota(snghbrs.begin(), snghbrs.end(), 0);
1043 // m_wallSendcounts: # cells to be sent to each rank
1044 m_wallSendcounts.resize(noRanks);
1045 // m_wallSendCells: local wallIds to be sent to each rank
1046 m_wallSendCells = maia::mpi::mpiExchangePointToPoint(&recvWallIds[0],
1047 &snghbrs[0],
1048 noRanks,
1049 &sendcounts[0],
1050 &snghbrs[0],
1051 noRanks,
1052 m_solver->m_StructuredComm,
1053 m_solver->domainId(),
1054 1,
1055 m_wallSendcounts.data());
1056
1057 // --> m_wallCellId2recvCell, m_wallSendCells and m_wallSendcounts can be used later to get information
1058 // from nearest wall neighbor, e.g., friction velocity at nearest wall cell
1059
1060 // Sanity checks
1061 // 1)
1062 if((signed)m_wallSendCells.size() != std::accumulate(&m_wallSendcounts[0], &m_wallSendcounts[0] + noRanks, 0))
1063 mTerm(1, "Neeeeeeiiiiiiiiin!");
1064
1065 // 2)
1066 if(cellmemSize != 0) {
1067 MInt min = std::numeric_limits<MInt>::max();
1068 MInt max = -1;
1069 for(auto id : m_wallSendCells) {
1070 // per wallDistInfo, we have one grid point more then grid cell
1071 if(id >= cellmemSize - noWallDistInfo) {
1072 cout << "DEBUG id=" << id << " cellmemSize=" << cellmemSize << endl;
1073 mTerm(1, "Something went wrong!");
1074 }
1075 min = mMin(id, min);
1076 max = mMax(id, max);
1077 }
1078 if(min != 0) {
1079 mTerm(1, "Scheisse: min!=0");
1080 }
1081 if(max != cellmemSize - 1 - noWallDistInfo) {
1082 cout << "SCHEISSE max=" << max << " cellmemSize=" << cellmemSize
1083 << " cellIdShifts[cellIdShifts.size()-2]=" << cellIdShifts[cellIdShifts.size() - 2] << endl;
1084 mTerm(1, "Scheisse ......");
1085 }
1086 }
1087
1088 // Some debug output
1089#ifndef NDEBUG
1090 cout << "::computeLocalWallDistance: rank=" << m_solver->domainId() << " cellmemSize=" << cellmemSize
1091 << " noTotalWallPoints=" << noTotalWallPoints << " #cells near wall=" << m_wallCellId2recvCell.size()
1092 << " #globalWallIds=" << globalWallIds.size() << endl;
1093#endif
1094}
maia::mpi::mpiExchangePointToPoint
std::vector< T > mpiExchangePointToPoint(const T *const sendBuffer, const MInt *const snghbrs, const MInt nosnghbrs, const MInt *const sendcounts, const MInt *const rnghbrs, const MInt nornghbrs, const MPI_Comm &mpi_comm, const MInt domainId, const MInt dataBlockSize, MInt *const recvcounts_=nullptr, MInt *const rdispls_=nullptr)
Definition: mpiexchange.h:1057

◆ computeWallDistances()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::computeWallDistances
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 329 of file fvstructuredbndrycnd2d.cpp.

329 {
330 vector<unique_ptr<StructuredWindowMap<nDim>>>& wallDistInfo = m_solver->m_windowInfo->m_wallDistInfoMap;
331 MInt noWallDistInfo = wallDistInfo.size(); // contains the size of the maps
332 MInt memSize = 0;
333
334 // initialize array with high numbers
335 for(MInt id = 0; id < m_noCells; ++id) {
336 m_cells->fq[FQ->WALLDISTANCE][id] = 99999;
337 }
338
339 // 1)
340 // memory will not be needed later ==> Scratch!
341
342 // 1.1) determine the storage size
343 // determine the size and to store the whole data
344 for(MInt i = 0; i < noWallDistInfo; ++i) {
345 MInt size = 1;
346 for(MInt dim = 0; dim < 2; ++dim) {
347 size *= (wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim] + 1);
348 }
349 memSize += size;
350 }
351
352 // 1.2) allocate the memory
353 MFloatScratchSpace coordMem(2 * memSize, AT_, "wallCoordinates");
354 MFloatPointerScratchSpace wallDistCoords(noWallDistInfo, AT_, "wallCoordsPointer");
355
356 MInt totMemSize = 0;
357 for(MInt i = 0; i < noWallDistInfo; ++i) {
358 MInt size = 1;
359 for(MInt dim = 0; dim < 2; ++dim) {
360 size *= (wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim] + 1);
361 }
362 wallDistCoords[i] = &coordMem[totMemSize];
363 totMemSize += 2 * size;
364 }
365
366 // 2)
367 // we do not need the corner points but the centre coordinates of the face
368 // we need to allocate the memory for the faces (==cell size) for the sponge face
369 // ->determine the size and allocate the memory (again only scratch)
370
371 // 2.1) calculate the number of cells to store.
372 MInt cellmemSize = 0;
373 // determine the size and to store the whole data
374 for(MInt i = 0; i < noWallDistInfo; ++i) {
375 MInt size = 1;
376 for(MInt dim = 0; dim < 2; ++dim) {
377 if(wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim] == 0) continue;
378 size *= (wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim]);
379 }
380 cellmemSize += size;
381 }
382 // 2.2) allocate the space for all the coordinates in Scratch!
383 // memory will not be needed later!
384 MFloatScratchSpace coordCellMem(2 * cellmemSize, AT_, "wallDistCellCoordinates");
385 MFloatPointerScratchSpace wallDistSurfCoords(noWallDistInfo, AT_, "wallDistCellCoordPointer");
386
387 MInt totCellMemSize = 0;
388 for(MInt i = 0; i < noWallDistInfo; ++i) {
389 MInt size = 1;
390 for(MInt dim = 0; dim < 2; ++dim) {
391 if(wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim] == 0) continue;
392 size *= (wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim]);
393 }
394 wallDistSurfCoords[i] = &coordCellMem[totCellMemSize];
395 totCellMemSize += 2 * size;
396 }
397
398 // 3) read in the coordinates of the grid points
399 // open file for reading the grid data
400 MString gridFileName = m_grid->m_gridInputFileName;
401
402 // open file and read number of solvers and uid
403 ParallelIoHdf5 pio(gridFileName, maia::parallel_io::PIO_READ, m_solver->m_StructuredComm);
404
405 // read the data in and distribute ist
406
407 // the split of reading and distributing is done on purpose!
408 // this is because we then know what the library is doing
409 // and not what the io_library such as hdf or netcdf should do
410 // but does not
411 // a good library should handle that directly! TEST IT
412 memSize = 1;
413 // read in the data if processor zero else read nothing!!!
414 ParallelIo::size_type ioOffset[2] = {0, 0};
415 ParallelIo::size_type ioSize[2] = {0, 0};
416 if(m_solver->domainId() == 0) {
417 for(MInt i = 0; i < noWallDistInfo; ++i) {
418 memSize = 1;
419 for(MInt dim = 1; dim >= 0; --dim) {
420 ioSize[dim] = (wallDistInfo[i]->end1[1 - dim] - wallDistInfo[i]->start1[1 - dim] + 1);
421 memSize *= ioSize[dim];
422 ioOffset[dim] = wallDistInfo[i]->start1[1 - dim];
423 }
424 // read in the data if processor zero else read nothing!!!
425 // determine the Solver name
426 MString bName = "block";
427 stringstream number;
428 number << wallDistInfo[i]->Id1;
429 bName += number.str();
430 pio.readArray(&wallDistCoords[i][0], bName, "x", 2, ioOffset, ioSize);
431 pio.readArray(&wallDistCoords[i][memSize], bName, "y", 2, ioOffset, ioSize);
432 }
433 } else {
434 for(MInt i = 0; i < noWallDistInfo; ++i) {
435 MString bName = "block";
436 stringstream number;
437 number << wallDistInfo[i]->Id1;
438 bName += number.str();
439 MFloat empty = 0;
440 pio.readArray(&empty, bName, "x", 2, ioOffset, ioSize);
441 pio.readArray(&empty, bName, "y", 2, ioOffset, ioSize);
442 }
443 }
444
445 // cout << "broadcasting wall-bc information" << endl;
446 // now broadcast the information to everyone!!!
447 MPI_Bcast(&wallDistCoords[0][0], totMemSize, MPI_DOUBLE, 0, m_solver->m_StructuredComm, AT_, "wallDistCoords[0][0]");
448 // cout << "broadcast end" << endl;
449
450 // 4) computing the coordinates of surface center from corner points;
451
452 for(MInt ii = 0; ii < noWallDistInfo; ++ii) {
453 MInt label, size1, count = 0;
454 for(label = 0; label < 2; label++) { // 2== dimensions
455 if(wallDistInfo[ii]->end1[label] - wallDistInfo[ii]->start1[label] == 0) break;
456 }
457 switch(label) {
458 case 0: {
459 size1 = wallDistInfo[ii]->end1[1] - wallDistInfo[ii]->start1[1] + 1;
460 for(MInt i = 0; i < size1 - 1; i++) {
461 MInt I = i;
462 MInt IP = i + 1;
463 wallDistSurfCoords[ii][count] = 0.5 * (wallDistCoords[ii][I] + wallDistCoords[ii][IP]);
464 wallDistSurfCoords[ii][count + (size1 - 1)] =
465 0.5 * (wallDistCoords[ii][I + size1] + wallDistCoords[ii][IP + size1]);
466 count++;
467 }
468 break;
469 }
470 case 1: {
471 size1 = wallDistInfo[ii]->end1[0] - wallDistInfo[ii]->start1[0] + 1;
472 for(MInt i = 0; i < size1 - 1; i++) {
473 MInt I = i;
474 MInt IP = i + 1;
475 wallDistSurfCoords[ii][count] = 0.5 * (wallDistCoords[ii][I] + wallDistCoords[ii][IP]);
476 wallDistSurfCoords[ii][count + (size1 - 1)] =
477 0.5 * (wallDistCoords[ii][I + size1] + wallDistCoords[ii][IP + size1]);
478 count++;
479 }
480 break;
481 }
482 default:
483 mTerm(1, AT_, "wall direction is messed up");
484 }
485 }
486
487 // now everyone can determine the shortest distance
488
489 m_log << "wall distance computation: searching for the nearest wall" << endl;
490 // cout << "seraching for the nearest" << endl;
491 // build a k-d-tree for a quick search:
492 // 1) rearragne the coordinates into points;
493 for(MInt i = 0; i < noWallDistInfo; ++i) {
494 MInt noPoints = 1;
495 for(MInt dim = 0; dim < 2; ++dim) {
496 if(wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim] == 0) continue;
497 noPoints *= (wallDistInfo[i]->end1[dim] - wallDistInfo[i]->start1[dim]);
498 }
499 // build up the points (surface centres)
500 vector<Point<2>> pts;
501 for(MInt j = 0; j < noPoints; ++j) {
502 Point<2> a(wallDistSurfCoords[i][j], wallDistSurfCoords[i][j + noPoints]);
503 pts.push_back(a);
504 }
505 // build up the tree
506 KDtree<2> tree(pts);
507 MFloat distance = -1.0;
508
509 // go through all the cells an determine the closest distance
510 for(MInt id = 0; id < m_noCells; ++id) {
511 distance = -1.1111111111111111; // to check
512 Point<2> pt(m_cells->coordinates[0][id], m_cells->coordinates[1][id]);
513 (void)tree.nearest(pt, distance);
514
515 // take minimum because another bc1000 might be further away than the current but would overwrite the actually
516 // closer distance
517 m_cells->fq[FQ->WALLDISTANCE][id] = mMin(m_cells->fq[FQ->WALLDISTANCE][id], distance);
518 }
519 }
520
521 m_log << "Wall Distance Computation SUCESSFUL: Saved minimum distance to next wall for all cells " << endl;
522}
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
maia::parallel_io::PIO_READ
const MInt PIO_READ
Definition: parallelio.h:40

◆ correctBndryCndIndices()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::correctBndryCndIndices
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2489 of file fvstructuredbndrycnd2d.cpp.

2489 {
2490 TRACE();
2491
2492 // in correcting cell Information
2493 for(MInt bcId = 0; bcId < m_noBndryCndIds; bcId++) {
2494 (this->*initBndryCndHandler[bcId])(bcId);
2495 }
2496}
StructuredBndryCnd< 2 >::initBndryCndHandler
BndryCndHandler * initBndryCndHandler
Definition: fvstructuredbndrycnd.h:191

◆ correctWallDistanceAtBoundary()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::correctWallDistanceAtBoundary ( MInt  bcId)

Definition at line 2515 of file fvstructuredbndrycnd2d.cpp.

2515 {
2516 // Besides the wall distance the ghost cells eventually need the utau of the nearest wall
2517 // neighbors; Here we assume that the correct nearest neighbor is chosen automatically
2518 // and that hence we don't need to apply any correction; we also don't apply any extrapolation
2519 // to the utau as is done with the wall distance
2520 if(std::find(FQ->fqNames.begin(), FQ->fqNames.end(), "wallDistance") != FQ->fqNames.end()) {
2521 const MInt IJ[nDim] = {1, m_nCells[1]};
2522 MInt* start = m_physicalBCMap[bcId]->start1;
2523 MInt* end = m_physicalBCMap[bcId]->end1;
2524
2525 // Here we find out the normal direction of the
2526 // boundary and the tangential direction.
2527 // This way we can make a general formulation of
2528 // the boundary condition
2529 const MInt face = m_physicalBCMap[bcId]->face;
2530 const MInt normalDir = face / 2;
2531 const MInt firstTangentialDir = (normalDir + 1) % nDim;
2532 const MInt normalDirStart = start[normalDir];
2533 const MInt firstTangentialStart = start[firstTangentialDir];
2534 const MInt firstTangentialEnd = end[firstTangentialDir];
2535 const MInt inc[nDim] = {IJ[normalDir], IJ[firstTangentialDir]};
2536 // determine indices for direction help
2537 const MInt n = (face % 2) * 2 - 1; //-1,+1
2538 const MInt g1 = normalDirStart + (MInt)(0.5 - (0.5 * (MFloat)n)); //+1,0
2539 const MInt g2 = normalDirStart + (MInt)(0.5 + (0.5 * (MFloat)n)); // 0,+1
2540 const MInt a1 = normalDirStart + (MInt)(0.5 - (1.5 * (MFloat)n)); //+2,-1
2541 const MInt a2 = normalDirStart + (MInt)(0.5 - (2.5 * (MFloat)n)); //+3,-2
2542
2543 for(MInt t1 = firstTangentialStart; t1 < firstTangentialEnd; t1++) {
2544 const MInt cellIdG1 = g1 * inc[0] + t1 * inc[1];
2545 const MInt cellIdG2 = g2 * inc[0] + t1 * inc[1];
2546 const MInt cellIdA1 = a1 * inc[0] + t1 * inc[1];
2547 const MInt cellIdA2 = a2 * inc[0] + t1 * inc[1];
2548
2549 m_cells->fq[FQ->WALLDISTANCE][cellIdG1] =
2550 F2 * m_cells->fq[FQ->WALLDISTANCE][cellIdA1] - m_cells->fq[FQ->WALLDISTANCE][cellIdA2];
2551 m_cells->fq[FQ->WALLDISTANCE][cellIdG2] =
2552 F2 * m_cells->fq[FQ->WALLDISTANCE][cellIdG1] - m_cells->fq[FQ->WALLDISTANCE][cellIdA1];
2553 }
2554 }
2555}

◆ distributeMapProperties()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::distributeMapProperties ( const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &  auxDataMap,
const std::vector< MInt > &  sendCells,
const std::vector< MInt > &  sendcounts,
const std::map< MInt, std::tuple< MInt, MInt, MFloat > > &  cellId2recvCell,
const std::vector< MInt > &  mapIndex,
MFloat * const  targetBuffer 
)
Date
2020-03-??

Definition at line 2049 of file fvstructuredbndrycnd2d.cpp.

2055 {
2056 static const MFloat UT = m_solver->m_Ma * sqrt(PV->TInfinity);
2057 MFloat* const RESTRICT p = &m_cells->pvariables[PV->P][0];
2058 MFloat* const RESTRICT rho = &m_cells->pvariables[PV->RHO][0];
2059
2060 // 1) Compute |cf| from cf-vector of al wall cells
2061 std::vector<MFloat> cf;
2062 MInt cnt = 0;
2063 for(MInt map = 0; map < (signed)auxDataMap.size(); ++map) {
2064 const MInt mapOffsetCf = m_cells->cfOffsets[mapIndex[map]];
2065 const MInt normalDir = auxDataMap[map]->face / 2;
2066 const MInt tangentialDir = 1 - normalDir;
2067 const MInt* const start = auxDataMap[map]->start1;
2068 const MInt* const end = auxDataMap[map]->end1;
2069 const MInt size = end[tangentialDir] - start[tangentialDir];
2070 ASSERT((signed)cf.size() == cnt, "cf.size() == cnt");
2071 cf.resize(cnt + size);
2072 MInt i = 0, j = 0;
2073 if((auxDataMap[map]->face) % 2 == 0) {
2074 i = start[normalDir]; // bottom wall
2075 } else {
2076 i = end[normalDir] - 1; // top wall //TODO_SS labels:FV,totest check later if the -1 is necessary
2077 }
2078 MInt& reali = (normalDir == 0) ? i : j;
2079 MInt& realj = (normalDir == 0) ? j : i;
2080 MInt ii = 0;
2081 for(j = start[tangentialDir]; j < end[tangentialDir]; ++j, ++cnt, ++ii) {
2082 // get id of boundary cell and cell above that one for extrapolation
2083 const MInt cellId = cellIndex(reali, realj);
2084 const MFloat T = m_solver->m_gamma * p[cellId] / rho[cellId];
2085 const MFloat nuLaminar = SUTHERLANDLAW(T) / rho[cellId];
2086 // for (MInt i = 0; i < size; ++i, ++cnt) {
2087 // The following is no longer cf
2088 cf[cnt] =
2089 sqrt(0.5
2090 * sqrt(POW2(m_cells->cf[mapOffsetCf + 0 * size + ii]) + POW2(m_cells->cf[mapOffsetCf + 1 * size + ii])))
2091 * UT / nuLaminar;
2092 }
2093 }
2094
2095 // 2) Fill sendBuffer
2096 ScratchSpace<MFloat> sendBuffer(std::max(1, (signed)sendCells.size()), FUN_, "sendBuffer");
2097 for(cnt = 0; cnt < (signed)sendCells.size(); ++cnt) {
2098 sendBuffer[cnt] = cf[sendCells[cnt]];
2099 }
2100
2101 // 3) Send & receive relevant data
2102 MInt noRanks;
2103 MPI_Comm_size(m_solver->m_StructuredComm, &noRanks);
2104 std::vector<MInt> snghbrs(noRanks);
2105 std::iota(snghbrs.begin(), snghbrs.end(), 0);
2106 std::vector<MFloat> recvBuffer = maia::mpi::mpiExchangePointToPoint(&sendBuffer[0],
2107 &snghbrs[0],
2108 noRanks,
2109 sendcounts.data(),
2110 &snghbrs[0],
2111 noRanks,
2112 m_solver->m_StructuredComm,
2113 m_solver->domainId(),
2114 1);
2115
2116 // 4) Scatter
2117 for(std::map<MInt, std::tuple<MInt, MInt, MFloat>>::const_iterator it = cellId2recvCell.begin();
2118 it != cellId2recvCell.end();
2119 ++it) {
2120 const MInt id1 = get<0>(it->second);
2121 const MInt id2 = get<1>(it->second);
2122 const MFloat weight = get<2>(it->second);
2123 const MFloat result = weight * recvBuffer[id1] + (1 - weight) * recvBuffer[id2];
2124 // TODO_SS labels:FV later maybe change this to yplus, also density correction misses
2125 targetBuffer[it->first] = result; // sqrt(0.5*result)*UT;//PV->UInfinity;
2126 }
2127}

◆ distributeWallAndFPProperties()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::distributeWallAndFPProperties
overridevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2013 of file fvstructuredbndrycnd2d.cpp.

2013 {
2014 // Store solid maps and fluid-porous maps in seperate lists
2015 vector<unique_ptr<StructuredWindowMap<nDim>>> wallDistInfo;
2016 vector<MInt> wallMapIndex;
2017 vector<unique_ptr<StructuredWindowMap<nDim>>> FPDistInfo;
2018 vector<MInt> FPMapIndex;
2019 MInt cnt = 0;
2020 for(auto& auxDataMap : m_auxDataMap) {
2021 unique_ptr<StructuredWindowMap<nDim>> temp = make_unique<StructuredWindowMap<nDim>>();
2022 m_solver->m_windowInfo->mapCpy(auxDataMap, temp);
2023 if(auxDataMap->isFluidPorousInterface) {
2024 FPDistInfo.push_back(move(temp));
2025 FPMapIndex.push_back(cnt);
2026 } else if((int)(auxDataMap->BC / 1000.0) == 1) {
2027 wallDistInfo.push_back(move(temp));
2028 wallMapIndex.push_back(cnt);
2029 }
2030 ++cnt;
2031 }
2032
2033 distributeMapProperties(wallDistInfo, m_wallSendCells, m_wallSendcounts, m_wallCellId2recvCell, wallMapIndex,
2034 m_cells->fq[FQ->UTAU]);
2035 distributeMapProperties(FPDistInfo, m_FPSendCells, m_FPSendcounts, m_FPCellId2recvCell, FPMapIndex,
2036 m_cells->fq[FQ->UTAU]);
2037 // Don't use FQ->UTAU2 -> use something which is only allocated in porous solver
2038 distributeMapProperties(FPDistInfo, m_FPSendCells_porous, m_FPSendcounts_porous, m_FPCellId2recvCell_porous,
2039 FPMapIndex, m_cells->fq[FQ->UTAU2]);
2040}
StructuredBndryCnd2D::distributeMapProperties
void distributeMapProperties(const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &, const std::vector< MInt > &, const std::vector< MInt > &, const std::map< MInt, std::tuple< MInt, MInt, MFloat > > &cellId2recvCell, const std::vector< MInt > &, MFloat *const)
Definition: fvstructuredbndrycnd2d.cpp:2049

◆ getCloserMap()

template<MBool isRans>
template<typename T >
void StructuredBndryCnd2D< isRans >::getCloserMap ( const MFloat * const  cells_fq_distance1,
std::vector< std::pair< MInt, MInt > > &  cellId2globalWallId1,
const MFloat * const  cells_fq_distance2,
std::vector< std::pair< MInt, MInt > > &  cellId2globalWallId2,
MFloat * const  cells_fq_distance,
T  comparator = {} 
)
Date
2020-03-19
Parameters
[in]cells_fq_distance1
[in/out]cellId2globalWallId1
[in]cells_fq_distance2
[in/out]cellId2globalWallId2
[out]cells_fq_distance
[in]comparator

Definition at line 1876 of file fvstructuredbndrycnd2d.cpp.

1881 {
1882 // Loop over all cells of current domain and decide which one is closer
1883 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1884 if(cellId2globalWallId1[cellId].first != -1 && cellId2globalWallId2[cellId].first != -1) {
1885 // In case both are possible candidates, apply comparator
1886 const MBool swtch = comparator(cellId, cells_fq_distance1[cellId], cells_fq_distance2[cellId]);
1887 if(swtch) {
1888 cellId2globalWallId2[cellId] = std::make_pair(-1, -1);
1889 cells_fq_distance[cellId] = cells_fq_distance1[cellId];
1890 } else {
1891 cellId2globalWallId1[cellId] = std::make_pair(-1, -1);
1892 cells_fq_distance[cellId] = cells_fq_distance2[cellId];
1893 }
1894 } else if(cellId2globalWallId1[cellId].first != -1) {
1895 cells_fq_distance[cellId] = cells_fq_distance1[cellId];
1896 } else {
1897 cells_fq_distance[cellId] = cells_fq_distance2[cellId];
1898 }
1899 }
1900}

◆ getPointIdFromCell()

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::getPointIdFromCell ( MInt  i,
MInt  j 
)
inline

Definition at line 2505 of file fvstructuredbndrycnd2d.cpp.

2505 {
2506 return i + (j * (m_nCells[1] + 1));
2507}

◆ getPointIdFromPoint()

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::getPointIdFromPoint ( MInt  origin,
MInt  incI,
MInt  incJ 
)
inline

Definition at line 2510 of file fvstructuredbndrycnd2d.cpp.

2510 {
2511 return origin + incI + incJ * m_nPoints[1];
2512}

◆ initBc1000()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc1000 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2558 of file fvstructuredbndrycnd2d.cpp.

2558 {
2559 (void)bcId;
2560}

◆ initBc1001()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc1001 ( MInt  )
inlinevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 60 of file fvstructuredbndrycnd2d.h.

60{};

◆ initBc1003()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc1003 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2563 of file fvstructuredbndrycnd2d.cpp.

2563 {
2564 (void)bcId;
2565 m_isothermalWallTemperature = F1;
2566 if(Context::propertyExists("isothermalWallTemperature", m_solverId)) {
2567 m_isothermalWallTemperature =
2568 Context::getSolverProperty<MFloat>("isothermalWallTemperature", m_solverId, AT_, &m_isothermalWallTemperature);
2569 }
2570}
Context::propertyExists
static MBool propertyExists(const MString &name, MInt solver=m_noSolvers)
This function checks if a property exists in general.
Definition: context.cpp:494

◆ initBc1004()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc1004 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2573 of file fvstructuredbndrycnd2d.cpp.

2573 { // moving adiabatic wall
2574 (void)bcId;
2575}

◆ initBc2001()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2001 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2578 of file fvstructuredbndrycnd2d.cpp.

2578 {
2579 correctWallDistanceAtBoundary(bcId);
2580}

◆ initBc2002()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2002 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2583 of file fvstructuredbndrycnd2d.cpp.

2583 {
2584 (void)bcId;
2585}

◆ initBc2003()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2003 ( MInt  )
inlinevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 65 of file fvstructuredbndrycnd2d.h.

65{};

◆ initBc2004()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2004 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2588 of file fvstructuredbndrycnd2d.cpp.

2588 {
2589 correctWallDistanceAtBoundary(bcId);
2590 bc2003(bcId);
2591}
StructuredBndryCnd2D::bc2003
void bc2003(MInt)
Subsonic in/outflow simple.
Definition: fvstructuredbndrycnd2d.cpp:3403

◆ initBc2005()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2005 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2594 of file fvstructuredbndrycnd2d.cpp.

2594 {
2595 (void)bcId;
2596}

◆ initBc2006()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2006 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2603 of file fvstructuredbndrycnd2d.cpp.

2603 {
2604 bc2003(bcId);
2605}

◆ initBc2007()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2007 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2609 of file fvstructuredbndrycnd2d.cpp.

2609 {
2610 correctWallDistanceAtBoundary(bcId);
2611}

◆ initBc2021()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2021 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2864 of file fvstructuredbndrycnd2d.cpp.

2864 {
2865 MInt* start = m_physicalBCMap[bcId]->start1;
2866 MInt* end = m_physicalBCMap[bcId]->end1;
2867 for(MInt var = 0; var < PV->noVariables; var++) {
2868 for(MInt i = start[0]; i < end[0]; i++) {
2869 for(MInt j = start[1]; j < end[1]; j++) {
2870 MInt cellId = cellIndex(i, j);
2871 MInt cellIdAdj = cellIndex(m_noGhostLayers, j);
2872 m_cells->pvariables[var][cellId] = m_cells->pvariables[var][cellIdAdj];
2873 }
2874 }
2875 }
2876 m_bc2021Gradient = 1.0;
2877 m_bc2021Gradient = Context::getSolverProperty<MFloat>("bc2021Gradient", m_solverId, AT_, &m_bc2021Gradient);
2878}

◆ initBc2199()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2199 ( MInt  )
inlinevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 71 of file fvstructuredbndrycnd2d.h.

71{};

◆ initBc2402()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2402 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2673 of file fvstructuredbndrycnd2d.cpp.

2673 {
2674 // for the channel flow we need the surface area of the entering and the outflow domain
2675 //==>determine Channelflow surface (assuming it does not change)
2676 MInt normal = 0;
2677 MInt Id = m_channelSurfaceIndexMap[bcId];
2678 if(Id < 0) cerr << "id smaller than zero ==> error" << endl;
2679 MInt* startface = &m_solver->m_windowInfo->channelSurfaceIndices[Id]->start1[0];
2680 MInt* endface = &m_solver->m_windowInfo->channelSurfaceIndices[Id]->end1[0];
2681 // find out which face it is
2682 MInt fixedind = -1;
2683
2684
2685 if(m_solver->m_initialCondition == 1233) {
2689 MFloat uTau = m_solver->m_ReTau * m_solver->m_Ma * sqrt(PV->TInfinity) / m_solver->m_Re;
2690 m_solver->m_deltaP = -CV->rhoInfinity * POW2(uTau) * F2 * (m_solver->m_channelLength) / m_solver->m_channelHeight;
2691 m_solver->m_channelPresInlet = PV->PInfinity;
2692 m_solver->m_channelPresOutlet = PV->PInfinity + m_solver->m_deltaP;
2693 m_log << "=========== Turb. Channel Flow BC Summary =========== " << endl;
2694 m_log << "-->Turbulent channel flow deltaP: " << m_solver->m_deltaP << endl;
2695 m_log << "-->channel friciton velocity: " << uTau << endl;
2696 m_log << "-->Channel pressure inflow: " << m_solver->m_channelPresInlet << endl;
2697 m_log << "-->Channel pressure outflow: " << m_solver->m_channelPresOutlet << endl;
2698 m_log << "=========== Turb. Channel Flow BC Summary Finished =========== " << endl;
2699 } else if(m_solver->m_initialCondition == 1234) {
2700 cout.precision(10);
2704
2705 m_solver->m_deltaP = -12.0 * PV->UInfinity * SUTHERLANDLAW(PV->TInfinity) * m_solver->m_channelLength
2706 / (POW2(m_solver->m_channelHeight) * m_solver->m_Re0);
2707 m_log << "=========== Lam. Channel Flow BC Summary =========== " << endl;
2708 m_log << "Laminar channel deltaP: " << m_solver->m_deltaP << endl;
2709 m_log << "Theoretical cD total (both channel walls): "
2710 << m_solver->m_deltaP * m_solver->m_channelHeight * m_solver->m_channelWidth
2711 / (0.5 * CV->rhoInfinity * POW2(PV->UInfinity))
2712 << endl;
2713 m_log << "Theoretical cD (single wall): "
2714 << m_solver->m_deltaP * m_solver->m_channelHeight * m_solver->m_channelWidth
2715 / (0.5 * CV->rhoInfinity * POW2(PV->UInfinity) * 2.0)
2716 << endl;
2717 m_log << "Theoretical cF total (both channel walls): "
2718 << m_solver->m_deltaP * m_solver->m_channelHeight
2719 / (0.5 * CV->rhoInfinity * POW2(PV->UInfinity) * m_solver->m_channelLength)
2720 << endl;
2721 m_log << "Theoretical cF (single wall): "
2722 << m_solver->m_deltaP * m_solver->m_channelHeight
2723 / (0.5 * CV->rhoInfinity * POW2(PV->UInfinity) * m_solver->m_channelLength * 2.0)
2724 << endl;
2725 m_log << "=========== Lam. Channel Flow BC Summary Finished =========== " << endl;
2726 }
2727
2728 for(MInt dim = 0; dim < nDim; dim++) {
2729 if(startface[dim] == endface[dim]) {
2730 // this is the normal
2731 fixedind = dim;
2732 if(startface[dim] == m_noGhostLayers) {
2733 normal = -1;
2734 } else {
2735 normal = 1;
2736 }
2737 }
2738 }
2739
2740 if(m_physicalBCMap[bcId]->face == 0) {
2741 MPI_Comm_rank(m_solver->m_commChannelIn, &m_channelInflowRank);
2742 }
2743
2744 switch(fixedind) {
2745 case 0: {
2746 MFloat surface = F0;
2747 MInt ii = 0;
2748 if(normal < 0) {
2749 ii = startface[0];
2750
2751 MInt cellId = 0;
2752 for(MInt j = startface[1]; j < endface[1]; j++) {
2753 cellId = cellIndex(ii, j);
2754 surface += sqrt(POW2(m_cells->cellMetrics[0][cellId]) + POW2(m_cells->cellMetrics[1][cellId]));
2755 }
2756 } else {
2757 ii = endface[0];
2758 // activate this or the method below with the 4 points for the surface calculation
2759 // this is less exact (for straight surf.) but works better for curved surfaces
2760
2761 MInt cellId = 0;
2762 for(MInt j = startface[1]; j < endface[1]; j++) {
2763 cellId = cellIndex(ii - 1, j);
2764 surface += sqrt(POW2(m_cells->cellMetrics[0][cellId]) + POW2(m_cells->cellMetrics[1][cellId]));
2765 }
2766 }
2767
2768 if(normal < 0) {
2769 MPI_Allreduce(&surface, &m_channelSurfaceIn, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelIn, AT_, "surface",
2770 "m_channelSurfaceIn");
2771 cout << "ChannelInSurface: " << m_channelSurfaceIn << endl;
2772 } else {
2773 MPI_Allreduce(&surface, &m_channelSurfaceOut, 1, MPI_DOUBLE, MPI_SUM, m_solver->m_commChannelOut, AT_,
2774 "surface", "m_channelSurfaceOut");
2775 cout << "ChannelOutSurface: " << m_channelSurfaceOut << endl;
2776 }
2777
2778 break;
2779 }
2780 case 1: {
2781 mTerm(1, AT_, "surface calculation for j faces(channel not implemented)");
2782 break;
2783 }
2784 default: {
2785 mTerm(1, AT_, "surface calculation for given faces(channel not implemented)");
2786 }
2787 }
2788}
StructuredCell::cellMetrics
MFloat ** cellMetrics
Definition: fvstructuredcell.h:42

◆ initBc2500()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2500 ( MInt  )
inlinevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 73 of file fvstructuredbndrycnd2d.h.

73{};

◆ initBc2501()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2501 ( MInt  )
inlinevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 74 of file fvstructuredbndrycnd2d.h.

74{};

◆ initBc2510()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2510 ( MInt  bcId)
virtual

Put values from field to gc to avoid nan, only matters for first Runge-Kutta Step at first timeStep

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2798 of file fvstructuredbndrycnd2d.cpp.

2798 {
2799 MInt* start = m_physicalBCMap[bcId]->start1;
2800 MInt* end = m_physicalBCMap[bcId]->end1;
2801
2802 for(MInt var = 0; var < PV->noVariables; var++) {
2803 for(MInt i = start[0]; i < end[0]; i++) {
2804 for(MInt j = start[1]; j < end[1]; j++) {
2805 MInt cellId = cellIndex(i, j);
2806 MInt cellIdAdj = cellIndex(m_noGhostLayers, j);
2807 m_cells->pvariables[var][cellId] = m_cells->pvariables[var][cellIdAdj];
2808 }
2809 }
2810 }
2811}

◆ initBc2511()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2511 ( MInt  )
inline

Definition at line 76 of file fvstructuredbndrycnd2d.h.

76{};

◆ initBc2600()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc2600 ( MInt  bcId)
virtual

In case of the initialStartup2600 flag is set (no matter if it was a restart or a initial startup) we load the values from the field to the GC and save them in the restart field

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2822 of file fvstructuredbndrycnd2d.cpp.

2822 {
2823 MInt* start = m_physicalBCMap[bcId]->start1;
2824 MInt* end = m_physicalBCMap[bcId]->end1;
2825
2826 m_solver->m_bc2600 = true;
2827
2828 switch(m_physicalBCMap[bcId]->face) {
2829 case 0: {
2830 if(m_solver->m_bc2600InitialStartup) {
2831 // First copy values from the field into the ghostcells
2832 for(MInt i = start[0]; i < end[0]; i++) {
2833 for(MInt j = start[1]; j < end[1]; j++) {
2834 const MInt cellId = cellIndex(i, j);
2835 const MInt cellIdAdj = cellIndex(m_noGhostLayers, j);
2836 for(MInt var = 0; var < PV->noVariables; var++) {
2837 m_cells->pvariables[var][cellId] = m_cells->pvariables[var][cellIdAdj];
2838 }
2839 }
2840 }
2841 }
2842
2843 // Then copy the values from the ghostcells into restart field
2844 for(MInt i = start[0]; i < end[0]; i++) {
2845 for(MInt j = m_noGhostLayers; j < end[1] - m_noGhostLayers; j++) {
2846 const MInt cellId = cellIndex(i, j);
2847 const MInt cellIdBc = i + (j - m_noGhostLayers) * m_noGhostLayers;
2848 for(MInt var = 0; var < PV->noVariables; var++) {
2849 m_solver->m_bc2600Variables[var][cellIdBc] = m_cells->pvariables[var][cellId];
2850 }
2851 }
2852 }
2853
2854 break;
2855 }
2856 default: {
2857 mTerm(1, AT_, "Face not implemented");
2858 break;
2859 }
2860 }
2861}

◆ initBc3000()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc3000 ( MInt  bcId)
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 2881 of file fvstructuredbndrycnd2d.cpp.

2881 {
2882 (void)bcId;
2883}

◆ initBc6002()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::initBc6002 ( MInt  )
inlineoverridevirtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 79 of file fvstructuredbndrycnd2d.h.

79{}; // Fluid-porous interface

◆ modifyFPDistance()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::modifyFPDistance ( const std::vector< std::unique_ptr< StructuredWindowMap< nDim > > > &  FPDistInfo,
MFloat * const  FP_distance,
std::vector< std::pair< MInt, MInt > > &  cellId2globalFPId,
const std::vector< MFloat > &  wallWeighting 
)

[Mößner, Michael, and Rolf Radespiel. "Modelling of turbulent flow over porous media using a volume averaging approach and a Reynolds stress model." Computers & Fluids 108 (2015): 25-42.]

Definition at line 1244 of file fvstructuredbndrycnd2d.cpp.

1247 {
1248 if(!m_solver->m_porous) return;
1249
1250 MFloat maxWallDistance = numeric_limits<MFloat>::max();
1251 maxWallDistance = Context::getSolverProperty<MFloat>("maxWallDistance", m_solverId, AT_, &maxWallDistance);
1252
1253 // If this is a porous solver, then discard the distance to fluid-porous interface
1254 if(m_solver->m_blockType == "porous") {
1255 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1256 FP_distance[cellId] = maxWallDistance;
1257 cellId2globalFPId[cellId] = std::make_pair(-1, -1);
1258 }
1259 // return;
1260 }
1261
1262 const MFloat* const RESTRICT por = &m_cells->fq[FQ->POROSITY][0];
1263 const MFloat* const RESTRICT Da = &m_cells->fq[FQ->DARCY][0];
1264 const auto& c_wd = m_solver->m_c_wd;
1265
1266 const MInt noWallDistInfo = FPDistInfo.size(); // contains the size of the maps
1267 std::vector<MInt> wallNormalIndex(noWallDistInfo);
1268 std::vector<MInt> normalDir(noWallDistInfo);
1269 std::vector<std::array<MInt, 2 * nDim>> startEndTangential(noWallDistInfo);
1270 // Number of cells on current rank to store
1271 MInt cellmemSize = 0;
1272 for(MInt nn = 0; nn < noWallDistInfo; ++nn) {
1273 const MInt* const start = FPDistInfo[nn]->start1;
1274 const MInt* const end = FPDistInfo[nn]->end1;
1275 const MInt face = FPDistInfo[nn]->face;
1276 normalDir[nn] = face / 2;
1277 wallNormalIndex[nn] = start[normalDir[nn]] - (face % 2);
1278
1279 MInt size = 1;
1280 for(MInt dim = 0; dim < nDim - 1; ++dim) {
1281 const MInt tangentialDir = (normalDir[nn] + (dim + 1)) % nDim;
1282 startEndTangential[nn][dim * 2 + 0] = start[tangentialDir]; // + m_noGhostLayers;
1283 startEndTangential[nn][dim * 2 + 1] = end[tangentialDir]; // - m_noGhostLayers;
1284 size *= startEndTangential[nn][dim * 2 + 1] - startEndTangential[nn][dim * 2 + 0];
1285 }
1286 cellmemSize += size;
1287#ifndef NDEBUG
1288 cout << "Debug output: FPDistInfo=" << nn << ": face=" << face << "; normalDir=" << normalDir[nn]
1289 << "; wallNormalIndex=" << wallNormalIndex[nn] << "; startEndTangential=" << startEndTangential[nn][0] << "|"
1290 << startEndTangential[nn][1] << "normal start/end=" << start[normalDir[nn]] << "|" << end[normalDir[nn]]
1291 << endl;
1292#endif
1293 }
1294
1295 // 2.2) allocate the space in Scratch!
1296 // memory will not be needed later!
1297 MFloatScratchSpace correctionMem(std::max(1, cellmemSize), AT_, "correctionMem");
1298
1299 // 4) computing the correction terms
1300 MInt cnt = 0;
1301 for(MInt nn = 0; nn < noWallDistInfo; ++nn) {
1302 MInt tIdx;
1303 MInt& i = (normalDir[nn] == 0) ? wallNormalIndex[nn] : tIdx;
1304 MInt& j = (normalDir[nn] == 0) ? tIdx : wallNormalIndex[nn];
1305 for(tIdx = startEndTangential[nn][0]; tIdx < startEndTangential[nn][1]; ++tIdx) {
1306 const MInt cellId = cellIndex(i, j);
1307 correctionMem[cnt++] = c_wd * sqrt(Da[cellId] / por[cellId]);
1308 }
1309 }
1310
1311 // Broadcast all to all ranks
1312 MInt noRanks;
1313 MPI_Comm_size(m_solver->m_StructuredComm, &noRanks);
1314 // recvcounts: How many wall surface coordinates to receive from each domain
1315 std::vector<MInt> recvcounts(noRanks);
1316 MPI_Allgather(&cellmemSize, 1, MPI_INT, recvcounts.data(), 1, MPI_INT, m_solver->m_StructuredComm, AT_, "cellmemSize",
1317 "recvcounts");
1318 std::vector<MInt> displs(noRanks);
1319 for(MInt r = 0; r < noRanks - 1; ++r) {
1320 displs[r + 1] = displs[r] + recvcounts[r];
1321 }
1322 const MInt noTotalWallPoints = displs[noRanks - 1] + recvcounts[noRanks - 1];
1323 // wallDistSurfCoordsAll: contains the wall coordinates of all wall surfaces of all ranks
1324 std::vector<MFloat> correctionMemAll(noTotalWallPoints);
1325 MPI_Allgatherv(&correctionMem[0], cellmemSize, MPI_DOUBLE, correctionMemAll.data(), recvcounts.data(), displs.data(),
1326 MPI_DOUBLE, m_solver->m_StructuredComm, AT_, "correctionMem", "correctionMemAll");
1327
1328 // IMPORTANT: Here we assume for simplicity that by the correction the nearest fluid-porous point does not change;
1329 // This is true if the material parameters along the fluid-porous interface are constant
1330 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1331 const MInt globalFPId1 = cellId2globalFPId[cellId].first;
1332 // If globalFPId1==-1 then also globalFPId2==-1 and if globalFPId1!=-1 then also globalFPId2!=-1
1333 if(globalFPId1 > -1) {
1334 const MInt globalFPId2 = cellId2globalFPId[cellId].second;
1335 const MFloat my_FP_distance = FP_distance[cellId];
1336 const MFloat correction1 = correctionMemAll[globalFPId1];
1337 const MFloat correction2 = correctionMemAll[globalFPId2];
1338 const MFloat correction = wallWeighting[cellId] * correction1 + (1 - wallWeighting[cellId]) * correction2;
1339 const MFloat my_new_FP_distance = my_FP_distance + correction;
1340 if(my_new_FP_distance < maxWallDistance) {
1341 FP_distance[cellId] = my_new_FP_distance;
1342 } else {
1343 FP_distance[cellId] = maxWallDistance;
1344 cellId2globalFPId[cellId] = std::make_pair(-1, -1);
1345 }
1346 }
1347 }
1348}

◆ pressure()

template<MBool isRans>
MFloat StructuredBndryCnd2D< isRans >::pressure ( MInt  cellId)
inline

Definition at line 4818 of file fvstructuredbndrycnd2d.cpp.

4818 {
4819 return m_cells->pvariables[PV->P][cellId];
4820}

◆ readAndDistributeSpongeCoordinates()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::readAndDistributeSpongeCoordinates

Definition at line 56 of file fvstructuredbndrycnd2d.cpp.

56 {
57 TRACE();
58
59 vector<unique_ptr<StructuredWindowMap<nDim>>>& spongeInfo = m_solver->m_windowInfo->m_spongeInfoMap;
60 MInt noSpongeInfo = spongeInfo.size(); // contains the size of the maps
61 MInt memSize = 0;
62
63 // 1)
64 // allocate the space for all the coordinates of the sponge in Scratch!
65 // memory will not be needed later ==> Scratch!
66
67 // 1.1) determine the storage size
68 // determine the size and to store the whole data
69 for(MInt i = 0; i < noSpongeInfo; ++i) {
70 MInt size = 1;
71 for(MInt dim = 0; dim < nDim; ++dim) {
72 size *= (spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim] + 1);
73 }
74 memSize += size;
75 }
76
77 // 1.2) allocate the memory
78 MFloatScratchSpace coordMem(nDim * memSize, AT_, "spongeCoordinates");
79 MFloatPointerScratchSpace spongeCoords(noSpongeInfo, AT_, "spongeCoordPointer");
80
81 MInt totMemSize = 0;
82 for(MInt i = 0; i < noSpongeInfo; ++i) {
83 MInt size = 1;
84 for(MInt dim = 0; dim < nDim; ++dim) {
85 size *= (spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim] + 1);
86 }
87 spongeCoords[i] = &coordMem[totMemSize];
88 totMemSize += nDim * size;
89 }
90
91
92 // 2)
93 // we do not need the corner points but the centre coordinates of the face
94 // we need to allocate the memory for the faces (==cell size) for the sponge face
95 // ->determine the size and allocate the memory (again only scratch)
96
97 // 2.1) calculate the number of cells to store.
98 MInt cellmemSize = 0;
99 // determine the size and to store the whole data
100 for(MInt i = 0; i < noSpongeInfo; ++i) {
101 MInt size = 1;
102 for(MInt dim = 0; dim < nDim; ++dim) {
103 if(spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim] == 0) continue;
104 size *= (spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim]);
105 }
106 cellmemSize += size;
107 }
108 // 2.2) allocate the space for all the coordinates in Scratch!
109 // memory will not be needed later!
110 MFloatScratchSpace coordCellMem(nDim * cellmemSize, AT_, "spongeCellCoordinates");
111 MFloatPointerScratchSpace spongeSurfCoords(noSpongeInfo, AT_, "spongeCellCoordPointer");
112
113 MInt totCellMemSize = 0;
114 for(MInt i = 0; i < noSpongeInfo; ++i) {
115 MInt size = 1;
116 for(MInt dim = 0; dim < nDim; ++dim) {
117 if(spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim] == 0) continue;
118 size *= (spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim]);
119 }
120 spongeSurfCoords[i] = &coordCellMem[totCellMemSize];
121 totCellMemSize += nDim * size;
122 }
123
124
125 // 3) read in the coordinates of the grid points
126 // open file for reading the grid data
127 MString gridFileName = m_grid->m_gridInputFileName;
128
129 // unique identifier needed to associate grid and solution in case of restart
130 // const char* aUID= new char[18];
131
132
133 // open file and read number of solvers and uid
134 ParallelIoHdf5 pio(gridFileName, maia::parallel_io::PIO_READ, m_StructuredComm);
135
136 // read the data in and distribute ist
137
138 // the split of reading and distributing is done on purpose!
139 // this is because we then know what the library is doing
140 // and not what the io_library such as hdf or netcdf should do
141 // but does not
142 // a good library should handle that directly! TEST IT
143 memSize = 1;
144 // read in the data if processor zero else read nothing!!!
145 ParallelIo::size_type ioOffset[2] = {0, 0};
146 ParallelIo::size_type ioSize[2] = {0, 0};
147 if(m_solver->domainId() == 0) {
148 for(MInt i = 0; i < noSpongeInfo; ++i) {
149 memSize = 1;
150 for(MInt dim = 1; dim >= 0; --dim) {
151 ioSize[dim] = (spongeInfo[i]->end1[1 - dim] - spongeInfo[i]->start1[1 - dim] + 1);
152 memSize *= ioSize[dim];
153 ioOffset[dim] = spongeInfo[i]->start1[1 - dim];
154 }
155 // read in the data if processor zero else read nothing!!!
156 // determine the Solver name
157 MString bName = "block";
158 stringstream number;
159 number << spongeInfo[i]->Id1;
160 bName += number.str();
161 pio.readArray(&spongeCoords[i][0], bName, "x", 2, ioOffset, ioSize);
162 pio.readArray(&spongeCoords[i][memSize], bName, "y", 2, ioOffset, ioSize);
163 }
164 } else {
165 for(MInt i = 0; i < noSpongeInfo; ++i) {
166 MString bName = "block";
167 stringstream number;
168 number << spongeInfo[i]->Id1;
169 bName += number.str();
170 MFloat empty = 0;
171 pio.readArray(&empty, bName, "x", 2, ioOffset, ioSize);
172 pio.readArray(&empty, bName, "y", 2, ioOffset, ioSize);
173 }
174 }
175
176 // now broadcast the information to everyone!!!
177 MPI_Bcast(&spongeCoords[0][0], totMemSize, MPI_DOUBLE, 0, m_StructuredComm, AT_, "spongeCoords[0][0]");
178
179 // 4) computing the coordinates of surface center from corner points;
180
181 for(MInt ii = 0; ii < noSpongeInfo; ++ii) {
182 MInt label, size1, count = 0;
183 for(label = 0; label < nDim; label++) { // 2== dimensions
184 if(spongeInfo[ii]->end1[label] - spongeInfo[ii]->start1[label] == 0) break;
185 }
186 switch(label) {
187 case 0: {
188 size1 = spongeInfo[ii]->end1[1] - spongeInfo[ii]->start1[1] + 1;
189 for(MInt i = 0; i < size1 - 1; i++) {
190 MInt I = i;
191 MInt IP = i + 1;
192 spongeSurfCoords[ii][count] = 0.5 * (spongeCoords[ii][I] + spongeCoords[ii][IP]);
193 spongeSurfCoords[ii][count + (size1 - 1)] =
194 0.5 * (spongeCoords[ii][I + size1] + spongeCoords[ii][IP + size1]);
195 count++;
196 }
197 break;
198 }
199 case 1: {
200 size1 = spongeInfo[ii]->end1[0] - spongeInfo[ii]->start1[0] + 1;
201 for(MInt i = 0; i < size1 - 1; i++) {
202 MInt I = i;
203 MInt IP = i + 1;
204 spongeSurfCoords[ii][count] = 0.5 * (spongeCoords[ii][I] + spongeCoords[ii][IP]);
205 spongeSurfCoords[ii][count + (size1 - 1)] =
206 0.5 * (spongeCoords[ii][I + size1] + spongeCoords[ii][IP + size1]);
207 count++;
208 }
209 break;
210 }
211 default:
212 mTerm(1, AT_, "sponge direction is messed up");
213 }
214 }
215
216 // now everyone can determine the shortest distance
217
218 m_log << "sponge layer build: searching for the nearest points (building simgma sponge)" << endl;
219 // cout << "seraching for the nearest points(building sigma sponge)" << endl;
220 // build a k-d-tree for a quick search:
221 // 1) rearragne the coordinates into points;
222 for(MInt i = 0; i < noSpongeInfo; ++i) {
223 MInt noPoints = 1;
224 for(MInt dim = 0; dim < nDim; ++dim) {
225 if(spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim] == 0) continue;
226 noPoints *= (spongeInfo[i]->end1[dim] - spongeInfo[i]->start1[dim]);
227 }
228 // build up the points (surface centres)
229 vector<Point<2>> pts;
230 for(MInt j = 0; j < noPoints; ++j) {
231 Point<2> a(spongeSurfCoords[i][j], spongeSurfCoords[i][j + noPoints]);
232 pts.push_back(a);
233 }
234
235 // build up the tree
236 KDtree<2> tree(pts);
237 MFloat distance = -1.0;
238 // go through all the cells an determine the closest distance
239 MFloat spongeThickness = spongeInfo[i]->spongeThickness;
240 for(MInt id = 0; id < m_noCells; ++id) {
241 distance = -1.1111111111111111; // to check
242 Point<2> pt(m_cells->coordinates[0][id], m_cells->coordinates[1][id]);
243 (void)tree.nearest(pt, distance);
244 if(distance <= spongeThickness) {
245 MFloat spongeFactor =
246 spongeInfo[i]->sigma * pow((spongeThickness - distance) / spongeThickness, spongeInfo[i]->beta);
247 m_cells->fq[FQ->SPONGE_FACTOR][id] = mMax(m_cells->fq[FQ->SPONGE_FACTOR][id], spongeFactor);
248 }
249 }
250 }
251
252 m_log << "Spone layer SUCESSFUL: finished building sigma sponge " << endl;
253}

◆ setUpNearMapComm()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::setUpNearMapComm ( const std::vector< std::pair< MInt, MInt > > &  cellId2globalWallId,
const std::vector< MInt > &  displs,
const std::vector< MFloat > &  weighting,
std::map< MInt, std::tuple< MInt, MInt, MFloat > > &  r_cellId2recvCell,
std::vector< MInt > &  r_sendcounts,
std::vector< MInt > &  r_sendCells 
)
Date
2020-03-19
Parameters
[in]cellId2globalWallId
[in]displs
[out]r_cellId2recvCell
[out]r_sendcounts
[out]r_sendCells

Definition at line 1914 of file fvstructuredbndrycnd2d.cpp.

1919 {
1920 // Determine unqiue set of globalWallIds
1921 std::set<MInt> globalWallIds;
1922 MInt noNearMapCells = 0;
1923 for(MInt cellId = 0; cellId < m_noCells; ++cellId, ++noNearMapCells) {
1924 if(cellId2globalWallId[cellId].first != -1) {
1925 globalWallIds.insert(cellId2globalWallId[cellId].first);
1926 // if cellId2globalWallId[cellId].first!=-1, then also cellId2globalWallId[cellId].second!=-1
1927 globalWallIds.insert(cellId2globalWallId[cellId].second);
1928 }
1929 }
1930
1931 const MInt noRanks = displs.size();
1932 auto getDomainId = [&displs, noRanks](const MInt id, const MInt hint = 0) {
1933 for(MInt rank = hint; rank < (signed)displs.size() - 1; ++rank) {
1934 if(displs[rank + 1] > id) return rank;
1935 }
1936 return noRanks - 1;
1937 };
1938
1939 // globalWallId2localId: mapping of globalWallId in wallDistSurfCoordsAll to local wallId on current rank
1940 std::vector<MInt> globalWallId2localId((globalWallIds.empty()) ? 0 : *globalWallIds.rbegin() + 1);
1941 // recvWallIds: local wallIds from each rank, requested by current rank
1942 std::vector<MInt> recvWallIds(globalWallIds.size());
1943 std::vector<MInt> sendcounts(noRanks, 0);
1944 MInt domainId_temp = 0;
1945 MInt localWallId = 0;
1946 for(auto it = globalWallIds.begin(); it != globalWallIds.end(); ++it, ++localWallId) {
1947 globalWallId2localId[*it] = localWallId;
1948 domainId_temp = getDomainId(*it, domainId_temp);
1949 // transform globalWallId "*it" to local wallId of rank domainId_temp
1950 recvWallIds[localWallId] = *it - displs[domainId_temp];
1951 ++sendcounts[domainId_temp];
1952 }
1953
1954 // r_cellId2recvCell: mapping from local cellId to the position, where its closest wall surface information is stored
1955 for(MInt cellId = 0; cellId < m_noCells; ++cellId) {
1956 // if cellId2globalWallId[cellId].first==-1, then also cellId2globalWallId[cellId].second==-1
1957 if(cellId2globalWallId[cellId].first == -1) continue;
1958 auto temp2 = std::make_tuple(globalWallId2localId[cellId2globalWallId[cellId].first],
1959 globalWallId2localId[cellId2globalWallId[cellId].second],
1960 weighting[cellId]);
1961 r_cellId2recvCell.insert({cellId, temp2});
1962 }
1963
1964 // Broadcast send/receive infos
1965 std::vector<MInt> snghbrs(noRanks);
1966 std::iota(snghbrs.begin(), snghbrs.end(), 0);
1967 // r_sendcounts: # cells to be sent to each rank
1968 r_sendcounts.resize(noRanks);
1969 // r_sendCells: local wallIds to be sent to each rank
1970 r_sendCells = maia::mpi::mpiExchangePointToPoint(&recvWallIds[0],
1971 &snghbrs[0],
1972 noRanks,
1973 &sendcounts[0],
1974 &snghbrs[0],
1975 noRanks,
1976 m_solver->m_StructuredComm,
1977 m_solver->domainId(),
1978 1,
1979 r_sendcounts.data());
1980
1981#ifndef NDEBUG
1982 // Sanity checks
1983 // 1)
1984 if((signed)r_sendCells.size() != std::accumulate(&r_sendcounts[0], &r_sendcounts[0] + noRanks, 0))
1985 mTerm(1, "Neeeeeeiiiiiiiiin!");
1986
1987 // 2) Check is not done for last rank
1988 const MInt cellmemSize =
1989 (m_solver->domainId() == noRanks - 1) ? 0 : displs[m_solver->domainId() + 1] - displs[m_solver->domainId()];
1990 if(cellmemSize != 0) {
1991 MInt min = std::numeric_limits<MInt>::max();
1992 MInt max = -1;
1993 for(auto id : r_sendCells) {
1994 if(id >= cellmemSize) mTerm(1, "Something went wrong!");
1995 min = mMin(id, min);
1996 max = mMax(id, max);
1997 }
1998 if(min != 0) {
1999 mTerm(1, "Scheisse: min!=0");
2000 }
2001 if(max != cellmemSize - 1) mTerm(1, "Scheisse......");
2002 }
2003
2004 // Some debug output
2005 cout << "::computeLocalWallDistance: rank=" << m_solver->domainId()
2006 << " cellmemSize=" << cellmemSize /*<< " noTotalWallPoints=" << noTotalWallPoints*/ << " #cells near wall="
2007 << r_cellId2recvCell.size() << " #globalWallIds=" << globalWallIds.size() << endl;
2008#endif
2009}

◆ shortestDistanceToLineElement()

template<MBool isRans>
MFloat StructuredBndryCnd2D< isRans >::shortestDistanceToLineElement ( const  MFloat(&)[nDim],
const  MFloat(&)[nDim],
const  MFloat(&)[nDim],
MFloat &  distance,
MFloat &  fraction 
)

Definition at line 1100 of file fvstructuredbndrycnd2d.cpp.

1104 {
1105 TRACE();
1106
1107 const MFloat linep1p2[nDim] = {p2[0] - p1[0], p2[1] - p1[1]};
1108 const MFloat linep1p2LengthPOW2 = POW2(linep1p2[0]) + POW2(linep1p2[1]);
1109 const MFloat linep1pT[nDim] = {pTarget[0] - p1[0], pTarget[1] - p1[1]};
1110 const MFloat dotProd = std::inner_product(std::begin(linep1p2), std::end(linep1p2), std::begin(linep1pT), 0.0);
1111 // If pProjection is required uncomment those lines
1112 // MFloat pProjection[nDim];
1113 if(dotProd < 0) {
1114 // std::copy(std::begin(p1), std::end(p1), std::begin(pProjection));
1115 distance = sqrt(std::inner_product(std::begin(linep1pT), std::end(linep1pT), std::begin(linep1pT), 0.0));
1116 fraction = 0.0;
1117 } else if(dotProd > linep1p2LengthPOW2) {
1118 // std::copy(std::begin(p2), std::end(p2), std::begin(pProjection));
1119 const MFloat linep2pT[nDim] = {pTarget[0] - p2[0], pTarget[1] - p2[1]};
1120 distance = sqrt(std::inner_product(std::begin(linep2pT), std::end(linep2pT), std::begin(linep2pT), 0.0));
1121 fraction = 1.0;
1122 } else {
1123 fraction = dotProd / linep1p2LengthPOW2;
1124 const MFloat pProjection[nDim] = {p1[0] + fraction * linep1p2[0], p1[1] + fraction * linep1p2[1]};
1125 const MFloat linepPpT[nDim] = {pTarget[0] - pProjection[0], pTarget[1] - pProjection[1]};
1126 distance = sqrt(std::inner_product(std::begin(linepPpT), std::end(linepPpT), std::begin(linepPpT), 0.0));
1127 }
1128 // fraction >0.5 would indicate that p2 is closer to pTarget then p1
1129 if(fraction < 0.0 || fraction > 1.0) mTerm(1, "fraction < 0.0 or fraction > 1.0");
1130 return linep1p2LengthPOW2;
1131}

◆ temperature()

template<MBool isRans>
MFloat StructuredBndryCnd2D< isRans >::temperature ( MInt  cellId)
inline

Definition at line 4823 of file fvstructuredbndrycnd2d.cpp.

4823 {
4824 const MFloat gamma = m_solver->m_gamma;
4825 MFloat t = gamma * m_cells->pvariables[PV->P][cellId] / m_cells->pvariables[PV->RHO][cellId];
4826 return t;
4827}

◆ updateSpongeLayer()

template<MBool isRans>
void StructuredBndryCnd2D< isRans >::updateSpongeLayer
virtual

Reimplemented from StructuredBndryCnd< 2 >.

Definition at line 256 of file fvstructuredbndrycnd2d.cpp.

256 {
257 // damp to infinity values
258 // for the moment only rho and rhoE are damped
259 const MFloat gammaMinusOne = m_solver->m_gamma - 1.0;
260 switch(m_spongeLayerType) {
261 case 1: {
262 MFloat deltaRhoE = F0, deltaRho = F0;
263 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
264 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
265 // compute the forcing term
266 // for the pressure or engery and the velocity
267 MInt cellId = cellIndex(i, j);
268 const MFloat rhoE =
269 m_cells->pvariables[PV->P][cellId] / gammaMinusOne
270 + F1B2 * m_cells->pvariables[PV->RHO][cellId]
271 * (POW2(m_cells->pvariables[PV->U][cellId]) + POW2(m_cells->pvariables[PV->V][cellId]));
272
273 deltaRhoE = rhoE - CV->rhoEInfinity;
274 deltaRho = m_cells->pvariables[PV->RHO][cellId] - (CV->rhoInfinity * m_targetDensityFactor);
275
276 m_cells->rightHandSide[CV->RHO_E][cellId] -= m_cells->fq[FQ->SPONGE_FACTOR][cellId] * deltaRhoE
277 * m_cells->cellJac[cellId]; // deltaP * m_cells->cellJac[cellId];
278 m_cells->rightHandSide[CV->RHO][cellId] -=
279 m_cells->fq[FQ->SPONGE_FACTOR][cellId] * deltaRho * m_cells->cellJac[cellId];
280 }
281 }
282 break;
283 }
284 case 2: {
285 // damp to values in the FQ field (set at startup, from RANS etc.)
286 // damp to values in the FQ field (set at startup, from RANS etc.)
287 MFloat deltaRhoE = F0, deltaRho = F0;
288 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
289 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
290 MInt cellId = cellIndex(i, j);
291 const MFloat rhoE =
292 m_cells->pvariables[PV->P][cellId] / gammaMinusOne
293 + F1B2 * m_cells->pvariables[PV->RHO][cellId]
294 * (POW2(m_cells->pvariables[PV->U][cellId]) + POW2(m_cells->pvariables[PV->V][cellId]));
295
296 deltaRhoE = rhoE - m_cells->fq[FQ->SPONGE_RHO_E][cellId];
297 deltaRho =
298 m_cells->pvariables[PV->RHO][cellId] - (m_cells->fq[FQ->SPONGE_RHO][cellId] * m_targetDensityFactor);
299 m_cells->rightHandSide[CV->RHO_E][cellId] -=
300 m_cells->fq[FQ->SPONGE_FACTOR][cellId] * deltaRhoE * m_cells->cellJac[cellId];
301 m_cells->rightHandSide[CV->RHO][cellId] -=
302 m_cells->fq[FQ->SPONGE_FACTOR][cellId] * deltaRho * m_cells->cellJac[cellId];
303 }
304 }
305
306 break;
307 }
308 case 6: {
309 // damp to PInfinity
310 const MFloat FgammaMinusOne = F1 / (m_solver->m_gamma - F1);
311 for(MInt j = m_noGhostLayers; j < m_nCells[0] - m_noGhostLayers; j++) {
312 for(MInt i = m_noGhostLayers; i < m_nCells[1] - m_noGhostLayers; i++) {
313 const MInt cellId = cellIndex(i, j);
314 const MFloat deltaP = (m_cells->pvariables[PV->P][cellId] - PV->PInfinity) * FgammaMinusOne;
315 m_cells->rightHandSide[CV->RHO_E][cellId] -=
316 m_cells->fq[FQ->SPONGE_FACTOR][cellId] * deltaP * m_cells->cellJac[cellId];
317 }
318 }
319 break;
320 }
321 default:
322 mTerm(1, AT_, "Sponge type doesn't exist");
323 }
324}
StructuredCell::rightHandSide
MFloat ** rightHandSide
Definition: fvstructuredcell.h:22

Friends And Related Function Documentation

◆ FvStructuredSolver2D

template<MBool isRans>
friend class FvStructuredSolver2D
friend

Definition at line 23 of file fvstructuredbndrycnd2d.h.

Member Data Documentation

◆ m_bc2021Gradient

template<MBool isRans>
MFloat StructuredBndryCnd2D< isRans >::m_bc2021Gradient
protected

Definition at line 152 of file fvstructuredbndrycnd2d.h.

◆ m_channelInflowRank

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::m_channelInflowRank
protected

Definition at line 144 of file fvstructuredbndrycnd2d.h.

◆ m_endCommChannel

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::m_endCommChannel
protected

Definition at line 143 of file fvstructuredbndrycnd2d.h.

◆ m_endCommPeriodic

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::m_endCommPeriodic
protected

Definition at line 141 of file fvstructuredbndrycnd2d.h.

◆ m_isothermalWallTemperature

template<MBool isRans>
MFloat StructuredBndryCnd2D< isRans >::m_isothermalWallTemperature
protected

Definition at line 151 of file fvstructuredbndrycnd2d.h.

◆ m_noPeriodicConnections

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::m_noPeriodicConnections
protected

Definition at line 147 of file fvstructuredbndrycnd2d.h.

◆ m_rescalingBLT

template<MBool isRans>
MFloat StructuredBndryCnd2D< isRans >::m_rescalingBLT
protected

Definition at line 150 of file fvstructuredbndrycnd2d.h.

◆ m_solver

template<MBool isRans>
FvStructuredSolver2D* StructuredBndryCnd2D< isRans >::m_solver
protected

Definition at line 139 of file fvstructuredbndrycnd2d.h.

◆ m_startCommChannel

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::m_startCommChannel
protected

Definition at line 142 of file fvstructuredbndrycnd2d.h.

◆ m_startCommPeriodic

template<MBool isRans>
MInt StructuredBndryCnd2D< isRans >::m_startCommPeriodic
protected

Definition at line 140 of file fvstructuredbndrycnd2d.h.

◆ nDim

template<MBool isRans>
constexpr const MInt StructuredBndryCnd2D< isRans >::nDim = 2
staticconstexpr

Definition at line 31 of file fvstructuredbndrycnd2d.h.

The documentation for this class was generated from the following files: