Skip to content

Commit

Permalink
Merge branch 'feature.bgkmovingaverage' into 'master.dev'
Browse files Browse the repository at this point in the history 
[feature.bgkmovingaverage] BGK Moving Average

See merge request piclas/piclas!743
  • Loading branch information
@scopplestone
scopplestone committed Dec 9, 2022
2 parents b951a3c + 2290f69 commit 44b5b89
Show file tree
Hide file tree
Showing 7 changed files with 263 additions and 228 deletions.
13 changes: 6 additions & 7 deletions docs/documentation/userguide/features-and-models/Bhatnagar-Gross-Krook.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,16 +80,15 @@ relaxation factors are below unity. The `BGK_DSMC_Ratio` gives the percentage of
utilized. In a couple BGK-DSMC simulation this variable indicates the boundary between BGK and DSMC. However, a value below 1 can
occur for pure BGK simulations due to low particle numbers, when an element is skipped.

<!-- An option is available to utilize a moving average for the variables used in the calculation of the relaxation frequency:
An option is available to utilize a moving average for the variables used in the calculation of the relaxation frequency:

Particles-BGK-MovingAverage = T

The purpose is to increase the sample size for steady gas flows. An extension of this feature to account for unsteady flows is to
limit the moving average to the last $N$ number of iterations, where the first value of the array is deleted and the most current
value is added at the end of the list:
The purpose is to increase the sample size and reduce the noise for steady gas flows. For this, a factor
0 < Particles-BGK-MovingAverageFac < 1 must be defined with which the old $$M^n$$ and newly sampled moments $$M$$ are weighted
to define the moments for the next time step $$M^{n+1}$$:

Particles-BGK-MovingAverageLength = 100
$$ M^{n+1}=Particles-BGK-MovingAverageFac*M+(1-Particles-BGK-MovingAverageFac)*M^n$$

Although this feature was tested with a hypersonic flow around a $70^\circ$ blunted cone and a nozzle expansion, a clear advantage
could not be observed, however, it might reduce the statistical noise for other application cases. -->
Particles-BGK-MovingAverageFac = 0.01

5 changes: 4 additions & 1 deletion regressioncheck/WEK_BGKFlow/Flow_N2_70degCone/parameter.ini
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,5 +118,8 @@ Particles-BGK-DoCellAdaptation = T
Particles-BGK-MinPartsPerCell = 12
Particles-BGK-SplittingDens = 3.8E20
! Couped BGK-DSMC simulation
Particles-CoupledBGKDSMC = F,T
Particles-CoupledBGKDSMC = F,T,T
Particles-BGK-DSMC-SwitchDens = 5E20

Particles-BGK-MovingAverage = F,F,T
nocrosscombination:Particles-CoupledBGKDSMC,Particles-BGK-MovingAverage
270 changes: 95 additions & 175 deletions src/particles/bgk/bgk_adaptation.f90
View file

Large diffs are not rendered by default.

78 changes: 73 additions & 5 deletions src/particles/bgk/bgk_colloperator.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,7 @@ MODULE MOD_BGK_CollOperator

CONTAINS

SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume)
! SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingPara, CorrectStep)
SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume, AveragingValues)
!===================================================================================================================================
!> Subroutine for the cell-local BGK collision operator:
!> 1.) Moment calculation: Summing up the relative velocities and their squares
Expand All @@ -54,7 +53,7 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume)
USE MOD_Particle_Vars ,ONLY: PartState, Species, PartSpecies, nSpecies, usevMPF, VarTimeStep
USE MOD_DSMC_Vars ,ONLY: SpecDSMC, DSMC, PartStateIntEn, PolyatomMolDSMC, RadialWeighting, CollInf
USE MOD_TimeDisc_Vars ,ONLY: dt
USE MOD_BGK_Vars ,ONLY: SpecBGK, BGKDoVibRelaxation!, BGKMovingAverageLength
USE MOD_BGK_Vars ,ONLY: SpecBGK, BGKDoVibRelaxation, BGKMovingAverage
USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor, BGK_MeanRelaxFactorCounter, BGK_MaxRelaxFactor, BGK_MaxRotRelaxFactor
USE MOD_BGK_Vars ,ONLY: BGK_PrandtlNumber
USE MOD_part_tools ,ONLY: GetParticleWeight
Expand All @@ -68,8 +67,7 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume)
REAL, INTENT(IN) :: NodeVolume
INTEGER, INTENT(INOUT) :: nPart
INTEGER, INTENT(INOUT) :: iPartIndx_Node(:)
! REAL, INTENT(INOUT), OPTIONAL :: AveragingPara(5,BGKMovingAverageLength)
! INTEGER, INTENT(INOUT), OPTIONAL :: CorrectStep
REAL, INTENT(INOUT), OPTIONAL :: AveragingValues(:)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -130,6 +128,9 @@ SUBROUTINE BGK_CollisionOperator(iPartIndx_Node, nPart, NodeVolume)
dens = totalWeight * Species(1)%MacroParticleFactor / NodeVolume
END IF

IF (BGKMovingAverage) THEN
CALL DoAveraging(dens, u2, u0ij, u2i, CellTemp, AveragingValues)
END IF
! Calculation of the rotational and vibrational degrees of freedom for molecules
nXiVibDOF=0 ! Initialize
IF (nSpecies.EQ.1) THEN
Expand Down Expand Up @@ -430,6 +431,73 @@ SUBROUTINE CalcMoments(nPart, iPartIndx_Node, nSpec, vBulkAll, totalWeight, tota

END SUBROUTINE CalcMoments

SUBROUTINE DoAveraging(dens, u2, u0ij, u2i, CellTemp, AverageValues)
!===================================================================================================================================
!> Moment calculation: Summing up the relative velocities and their squares
!===================================================================================================================================
! MODULES
USE MOD_Particle_Vars ,ONLY: Species
USE MOD_BGK_Vars ,ONLY: BGKCollModel, BGKMovingAverageFac
USE MOD_Globals_Vars ,ONLY: BoltzmannConst
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
REAL, INTENT(INOUT) :: dens, u0ij(3,3), u2i(3), u2, CellTemp, AverageValues(:)
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES

!===================================================================================================================================
IF (BGKCollModel.EQ.1) THEN
IF (AverageValues(1).EQ.0.0) THEN
AverageValues(1) = dens
AverageValues(2) = u2
AverageValues(3) = u0ij(1,1); AverageValues(4) = u0ij(2,2); AverageValues(5) = u0ij(3,3);
AverageValues(6) = u0ij(1,2); AverageValues(7) = u0ij(1,3); AverageValues(8) = u0ij(2,3);
ELSE
AverageValues(1) = BGKMovingAverageFac*dens + (1.-BGKMovingAverageFac)*AverageValues(1)
AverageValues(2) = BGKMovingAverageFac*u2 + (1.-BGKMovingAverageFac)*AverageValues(2)
AverageValues(3) = BGKMovingAverageFac*u0ij(1,1) + (1.-BGKMovingAverageFac)*AverageValues(3)
AverageValues(4) = BGKMovingAverageFac*u0ij(2,2) + (1.-BGKMovingAverageFac)*AverageValues(4)
AverageValues(5) = BGKMovingAverageFac*u0ij(3,3) + (1.-BGKMovingAverageFac)*AverageValues(5)
AverageValues(6) = BGKMovingAverageFac*u0ij(1,2) + (1.-BGKMovingAverageFac)*AverageValues(6)
AverageValues(7) = BGKMovingAverageFac*u0ij(1,3) + (1.-BGKMovingAverageFac)*AverageValues(7)
AverageValues(8) = BGKMovingAverageFac*u0ij(2,3) + (1.-BGKMovingAverageFac)*AverageValues(8)
END IF
dens = AverageValues(1)
u2 = AverageValues(2)
u0ij(1,1)=AverageValues(3); u0ij(2,2)=AverageValues(4); u0ij(3,3)=AverageValues(5);
u0ij(1,2)=AverageValues(6); u0ij(1,3)=AverageValues(7); u0ij(2,3)=AverageValues(8);
ELSE IF (BGKCollModel.EQ.2) THEN
IF (AverageValues(1).EQ.0.0) THEN
AverageValues(1) = dens
AverageValues(2) = u2
AverageValues(3:5) = u2i(1:3)
ELSE
AverageValues(1) = BGKMovingAverageFac*dens + (1.-BGKMovingAverageFac)*AverageValues(1)
AverageValues(2) = BGKMovingAverageFac*u2 + (1.-BGKMovingAverageFac)*AverageValues(2)
AverageValues(3:5) = BGKMovingAverageFac*u2i(1:3) + (1.-BGKMovingAverageFac)*AverageValues(3:5)
END IF
dens = AverageValues(1)
u2 = AverageValues(2)
u2i(1:3) = AverageValues(3:5)
ELSE IF (BGKCollModel.EQ.3) THEN
IF (AverageValues(1).EQ.0.0) THEN
AverageValues(1) = dens
AverageValues(2) = u2
ELSE
AverageValues(1) = BGKMovingAverageFac*dens + (1.-BGKMovingAverageFac)*AverageValues(1)
AverageValues(2) = BGKMovingAverageFac*u2 + (1.-BGKMovingAverageFac)*AverageValues(2)
END IF
dens = AverageValues(1)
u2 = AverageValues(2)
END IF
CellTemp = Species(1)%MassIC * u2 / (3.0*BoltzmannConst)

END SUBROUTINE DoAveraging

SUBROUTINE CalcInnerDOFs(nSpec, EVibSpec, ERotSpec, totalWeightSpec, TVibSpec, TRotSpec, InnerDOF, Xi_VibSpec, Xi_Vib_oldSpec &
,Xi_RotSpec)
!===================================================================================================================================
Expand Down
88 changes: 72 additions & 16 deletions src/particles/bgk/bgk_init.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,9 +68,8 @@ SUBROUTINE DefineParametersBGK()
'cell refinement')
CALL prms%CreateLogicalOption('Particles-BGK-MovingAverage', 'Enable a moving average of variables for the calculation '//&
'of the cell temperature for relaxation frequencies','.FALSE.')
CALL prms%CreateIntOption( 'Particles-BGK-MovingAverageLength', 'Use the moving average with a fixed array length, where '//&
'the first values are dismissed and the last values updated '//&
'with current iteration to account for for unsteady flows','1')
CALL prms%CreateRealOption( 'Particles-BGK-MovingAverageFac', 'Use the moving average of moments M with '//&
'M^n+1=AverageFac*M+(1-AverageFac)*M^n','0.01')
CALL prms%CreateRealOption( 'Particles-BGK-SplittingDens', 'Octree-refinement will only be performed above this number '//&
'density', '0.0')
CALL prms%CreateLogicalOption('Particles-BGK-DoVibRelaxation', 'Enable modelling of vibrational excitation','.TRUE.')
Expand Down Expand Up @@ -165,12 +164,10 @@ SUBROUTINE InitBGK()
! Moving Average
BGKMovingAverage = GETLOGICAL('Particles-BGK-MovingAverage')
IF(BGKMovingAverage) THEN
CALL abort(__STAMP__,' ERROR BGK Init: Moving average is currently not implemented!')
BGKMovingAverageLength = GETINT('Particles-BGK-MovingAverageLength')
BGKMovingAverageFac= GETREAL('Particles-BGK-MovingAverageFac')
IF ((BGKMovingAverageFac.LE.0.0).OR.(BGKMovingAverageFac.GE.1.)) CALL abort(__STAMP__,'Particles-BGK-MovingAverageFac must be between 0 and 1!')
CALL BGK_init_MovingAverage()
IF(RadialWeighting%DoRadialWeighting.OR.VarTimeStep%UseVariableTimeStep) THEN
CALL abort(__STAMP__,' ERROR BGK Init: Moving average is neither implemented with radial weighting nor variable time step!')
END IF
IF(nSpecies.GT.1) CALL abort(__STAMP__,'nSpecies >1 and molecules not implemented for BGK averaging!')
END IF
IF(MoleculePresent) THEN
! Vibrational modelling
Expand Down Expand Up @@ -200,7 +197,7 @@ SUBROUTINE BGK_init_MovingAverage()
!> Initialization of the arrays for the sampling of the moving average
!===================================================================================================================================
! MODULES
USE MOD_BGK_Vars ,ONLY: ElemNodeAveraging, BGKMovingAverageLength
USE MOD_BGK_Vars ,ONLY: ElemNodeAveraging, BGKCollModel
USE MOD_Mesh_Vars ,ONLY: nElems
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
Expand All @@ -215,13 +212,14 @@ SUBROUTINE BGK_init_MovingAverage()
ALLOCATE(ElemNodeAveraging(nElems))
DO iElem = 1, nElems
ALLOCATE(ElemNodeAveraging(iElem)%Root)
IF (BGKMovingAverageLength.GT.1) THEN
ALLOCATE(ElemNodeAveraging(iElem)%Root%AverageValues(5,BGKMovingAverageLength))
ElemNodeAveraging(iElem)%Root%AverageValues = 0.0
ELSE
ALLOCATE(ElemNodeAveraging(iElem)%Root%AverageValues(5,1))
ElemNodeAveraging(iElem)%Root%AverageValues = 0.0
IF (BGKCollModel.EQ.1) THEN
ALLOCATE(ElemNodeAveraging(iElem)%Root%AverageValues(8))
ELSE IF (BGKCollModel.EQ.2) THEN
ALLOCATE(ElemNodeAveraging(iElem)%Root%AverageValues(5))
ELSE IF (BGKCollModel.EQ.3) THEN
ALLOCATE(ElemNodeAveraging(iElem)%Root%AverageValues(2))
END IF
ElemNodeAveraging(iElem)%Root%AverageValues = 0.0
ElemNodeAveraging(iElem)%Root%CorrectStep = 0
END DO

Expand All @@ -246,9 +244,67 @@ SUBROUTINE FinalizeBGK()
!===================================================================================================================================

SDEALLOCATE(SpecBGK)
SDEALLOCATE(ElemNodeAveraging)
SDEALLOCATE(BGK_QualityFacSamp)
IF(BGKMovingAverage) CALL DeleteElemNodeAverage()

END SUBROUTINE FinalizeBGK

SUBROUTINE DeleteElemNodeAverage()
!----------------------------------------------------------------------------------------------------------------------------------!
! Delete the pointer tree ElemNodeVol
!----------------------------------------------------------------------------------------------------------------------------------!
! MODULES !
!----------------------------------------------------------------------------------------------------------------------------------!
USE MOD_Globals
USE MOD_BGK_Vars
USE MOD_DSMC_Vars ,ONLY: DSMC
USE MOD_Mesh_Vars ,ONLY: nElems
!----------------------------------------------------------------------------------------------------------------------------------!
IMPLICIT NONE
! INPUT VARIABLES
!----------------------------------------------------------------------------------------------------------------------------------!
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER :: iElem
!===================================================================================================================================
IF(.NOT.DSMC%UseOctree) RETURN
DO iElem=1,nElems
CALL DeleteNodeAverage(ElemNodeAveraging(iElem)%Root)
DEALLOCATE(ElemNodeAveraging(iElem)%Root)
END DO
DEALLOCATE(ElemNodeAveraging)
END SUBROUTINE DeleteElemNodeAverage


RECURSIVE SUBROUTINE DeleteNodeAverage(NodeAverage)
!----------------------------------------------------------------------------------------------------------------------------------!
! Check if the Node has subnodes and delete them
!----------------------------------------------------------------------------------------------------------------------------------!
! MODULES !
!----------------------------------------------------------------------------------------------------------------------------------!
USE MOD_Globals
USE MOD_BGK_Vars
USE MOD_DSMC_Vars ,ONLY: DSMC
USE MOD_Particle_Vars ,ONLY: Symmetry
!----------------------------------------------------------------------------------------------------------------------------------!
IMPLICIT NONE
! INPUT VARIABLES
TYPE (tNodeAverage), INTENT(INOUT) :: NodeAverage
!----------------------------------------------------------------------------------------------------------------------------------!
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER :: iLoop, nLoop
!===================================================================================================================================
nLoop = 2**Symmetry%Order
IF(ASSOCIATED(NodeAverage%SubNode)) THEN
DO iLoop = 1, nLoop
CALL DeleteNodeAverage(NodeAverage%SubNode(iLoop))
END DO
DEALLOCATE(NodeAverage%SubNode)
END IF

END SUBROUTINE DeleteNodeAverage

END MODULE MOD_BGK_Init
24 changes: 10 additions & 14 deletions src/particles/bgk/bgk_main.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ SUBROUTINE BGK_DSMC_main(stage_opt)
USE MOD_BGK_Adaptation ,ONLY: BGK_octree_adapt, BGK_quadtree_adapt
USE MOD_Particle_Vars ,ONLY: PEM, Species, WriteMacroVolumeValues, Symmetry, usevMPF
USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation,BGKDSMCSwitchDens
! USE MOD_BGK_Vars ,ONLY: BGKMovingAverage,ElemNodeAveraging,BGKMovingAverageLength
USE MOD_BGK_Vars ,ONLY: BGKMovingAverage,ElemNodeAveraging
USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor,BGK_QualityFacSamp
USE MOD_BGK_Vars ,ONLY: BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber
USE MOD_BGK_CollOperator ,ONLY: BGK_CollisionOperator
Expand Down Expand Up @@ -130,13 +130,11 @@ SUBROUTINE BGK_DSMC_main(stage_opt)
BGK_MeanRelaxFactorCounter = 0; BGK_MeanRelaxFactor = 0.; BGK_MaxRelaxFactor = 0.; BGK_MaxRotRelaxFactor = 0.
BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0.
END IF
! IF (BGKMovingAverage) THEN
! CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID), &
! ElemNodeAveraging(iElem)%Root%AverageValues(1:5,1:BGKMovingAverageLength), &
! CorrectStep = ElemNodeAveraging(iElem)%Root%CorrectStep)
! ELSE
IF (BGKMovingAverage) THEN
CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID), ElemNodeAveraging(iElem)%Root%AverageValues(:))
ELSE
CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID))
! END IF
END IF
DEALLOCATE(iPartIndx_Node)
IF(DSMC%CalcQualityFactors) THEN
IF((Time.GE.(1-DSMC%TimeFracSamp)*TEnd).OR.WriteMacroVolumeValues) THEN
Expand Down Expand Up @@ -169,7 +167,7 @@ SUBROUTINE BGK_main(stage_opt)
USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem
USE MOD_BGK_Adaptation ,ONLY: BGK_octree_adapt, BGK_quadtree_adapt
USE MOD_Particle_Vars ,ONLY: PEM, WriteMacroVolumeValues, WriteMacroSurfaceValues, Symmetry, DoVirtualCellMerge, VirtMergedCells
USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation!, BGKMovingAverage, ElemNodeAveraging, BGKMovingAverageLength
USE MOD_BGK_Vars ,ONLY: DoBGKCellAdaptation, BGKMovingAverage, ElemNodeAveraging
USE MOD_BGK_Vars ,ONLY: BGK_MeanRelaxFactor,BGK_MeanRelaxFactorCounter,BGK_MaxRelaxFactor,BGK_QualityFacSamp
USE MOD_BGK_Vars ,ONLY: BGK_MaxRotRelaxFactor, BGK_PrandtlNumber, BGK_ExpectedPrandtlNumber
USE MOD_BGK_CollOperator ,ONLY: BGK_CollisionOperator
Expand Down Expand Up @@ -271,13 +269,11 @@ SUBROUTINE BGK_main(stage_opt)
BGK_PrandtlNumber=0.; BGK_ExpectedPrandtlNumber=0.
END IF

! IF (BGKMovingAverage) THEN
! CALL BGK_CollisionOperator(iPartIndx_Node, nPart, ElemVolume_Shared(CNElemID), &
! ElemNodeAveraging(iElem)%Root%AverageValues(1:5,1:BGKMovingAverageLength), &
! CorrectStep = ElemNodeAveraging(iElem)%Root%CorrectStep)
! ELSE
IF (BGKMovingAverage) THEN
CALL BGK_CollisionOperator(iPartIndx_Node, nPartMerged, elemVolume,ElemNodeAveraging(iElem)%Root%AverageValues(:))
ELSE
CALL BGK_CollisionOperator(iPartIndx_Node, nPartMerged, elemVolume)
! END IF
END IF
DEALLOCATE(iPartIndx_Node)
IF(DSMC%CalcQualityFactors) THEN
IF((Time.GE.(1-DSMC%TimeFracSamp)*TEnd).OR.WriteMacroVolumeValues) THEN
Expand Down
Loading

0 comments on commit 44b5b89

Please sign in to comment.