Wrong number density calculation in axisymmetric DSMC case

[MaXXe4th]

Hello everyone,

Case files: piclas_issue.zip

Output number density

I am simulating the exhaust plume of an electrical thruster in a vacuum chamber using PICLas. For this I want to first simulate a neutral Ar stream in DSMC which will be my background gas in future PIC simulations.

I use an axisymmetric domain (x_max 0.4 cm, y_max 0.4 cm) and used tutorial 7.6 as my template.

The domain has a residual pressure of 0.2 Pa which translates to about 5E19 $m^{-3}$ number density which I want to keep steady inside the domain during the simulation.

The left side (x=0) consists of wall boundary conditions and a small inlet area (the thruster). Bottom is symmetric_axis and front and back are symmetric_dim. Right now I have ungrouped the top and right boundaries, in order to experiment with the conditions. Essentially I want them both to be open condition, but keep the pressure stable at 0.2 Pa. (Adaptive boundaries)

The issue

Whatever I try, the number density at the top boundary condition decreases, while the right works perfectly fine. As I use the same conditions for the right outlet, my guess is that the big surface area of the top boundary (which is effectively an cylindric plane at radius 0.4) is causing this.

What I've tried

• If I set the top boundary as wall, everything works fine. However I would expect open conditions with steady pressure to also work. As this is closer to my application I would like to keep open boundaries.
• I tried a Surface Flux with fixed PartDensity 4.83E19 and 0.0 VeloIC at the top boundary
• I experimenterd with different radial weighting. In the current state I have multiple 100 particles in the topmost cells
• I set SurfFluxSides to different values but didn't see any difference in the behavior. Right now I have it at 200.

For the Adaptive pressure boundary I use

Part-Species1-Surfaceflux2-Adaptive = TRUE
Part-Species1-Surfaceflux2-Adaptive-Type = 1
Part-Species1-Surfaceflux2-Adaptive-Pressure = 0.2

and

AdaptiveBC-TruncateRunningAverage = T ! DEFAULT: F
AdaptiveBC-SamplingIteration = 100

I would be grateful for any insight into what the problem is. If you need further information, please let me know. (m.mudra@unibw.de)

Best regards

Max

[pnizenkov]

Hi Max,

it appears to be a bug that we recently fixed. We are currently preparing the next release 3.7.0, where the bugfix will be included. If you like, you can add this else case in src/particles/emission/particle_surface_flux.f90 below line 892:

  IF(Species(iSpec)%Surfaceflux(iSF)%nVFRSub(iSide,1).GT.0.) THEN
    ! case orthogonal to symmetry axis
  ELSE
    CALL RANDOM_NUMBER(RandVal)
    PartInsSubSide = INT(Species(iSpec)%Surfaceflux(iSF)%PartDensity / Species(iSpec)%MacroParticleFactor &
      * dtVar*RKdtFrac * Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%nVFR + RandVal)
  END IF

Cheers,
Paul

[MaXXe4th]

Hey Paul,

thanks for the quick response. I am using 3.6.0 and the source file seems to already have exactly this ELSE case. I compiled piclas from source, so I would think this file has been used.

particle_surface_flux.f90:

! Species-specific time step
IF(VarTimeStep%UseSpeciesSpecific) THEN
  dtVar = dt * Species(iSpec)%TimeStepFactor
ELSE
  dtVar = dt
END IF

IF(ParticleWeighting%UseSubdivision) THEN
  PartInsSubSide = 0
  IF(Species(iSpec)%Surfaceflux(iSF)%nVFRSub(iSide,1).GT.0.) THEN
    DO iSub = 1, ParticleWeighting%nSubSides
      CALL RANDOM_NUMBER(RandVal)
      ParticleWeighting%PartInsSide(iSub) = INT(Species(iSpec)%Surfaceflux(iSF)%PartDensity / Species(iSpec)%MacroParticleFactor &
              * dtVar*RKdtFrac * Species(iSpec)%Surfaceflux(iSF)%nVFRSub(iSide,iSub)+ RandVal)
      PartInsSubSide = PartInsSubSide + ParticleWeighting%PartInsSide(iSub)
    END DO
  END IF
ELSE
  CALL RANDOM_NUMBER(RandVal)
  PartInsSubSide = INT(Species(iSpec)%Surfaceflux(iSF)%PartDensity / Species(iSpec)%MacroParticleFactor &
    * dtVar*RKdtFrac * Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%nVFR + RandVal)
END IF

END SUBROUTINE CalcPartInsVarWeight

Best regards

Max

[pnizenkov]

Yeah, but in your case of an axisymmetric simulation, ParticleWeighting%UseSubdivision is enabled. And you need an ELSE to the IF(Species(iSpec)%Surfaceflux(iSF)%nVFRSub(iSide,1).GT.0.) THEN case.

[MaXXe4th]

A, I see. Maybe the numbering got moved around a bit. I'll try that out and let you know if that worked.

[MaXXe4th]

Hey Paul,

I made the changes to the file:

IF(ParticleWeighting%UseSubdivision) THEN
  PartInsSubSide = 0
  IF(Species(iSpec)%Surfaceflux(iSF)%nVFRSub(iSide,1).GT.0.) THEN
    DO iSub = 1, ParticleWeighting%nSubSides
      CALL RANDOM_NUMBER(RandVal)
      ParticleWeighting%PartInsSide(iSub) = INT(Species(iSpec)%Surfaceflux(iSF)%PartDensity / Species(iSpec)%MacroParticleFactor &
              * dtVar*RKdtFrac * Species(iSpec)%Surfaceflux(iSF)%nVFRSub(iSide,iSub)+ RandVal)
      PartInsSubSide = PartInsSubSide + ParticleWeighting%PartInsSide(iSub)
    END DO
  ELSE
    CALL RANDOM_NUMBER(RandVal)
    PartInsSubSide = INT(Species(iSpec)%Surfaceflux(iSF)%PartDensity / Species(iSpec)%MacroParticleFactor &
      * dtVar*RKdtFrac * Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%nVFR + RandVal)
  END IF
ELSE
  CALL RANDOM_NUMBER(RandVal)
  PartInsSubSide = INT(Species(iSpec)%Surfaceflux(iSF)%PartDensity / Species(iSpec)%MacroParticleFactor &
    * dtVar*RKdtFrac * Species(iSpec)%Surfaceflux(iSF)%SurfFluxSubSideData(iSample,jSample,iSide)%nVFR + RandVal)
END IF

I created a new build directory ./piclas/build_DSMC_2 and used these settings:

cmake .. -DPICLAS_TIMEDISCMETHOD=DSMC -DPICLAS_EQNSYSNAME=maxwell

I created symbolic links:

ln -s $PICLAS_PATH/build_DSMC_2/bin/piclas
ln -s $PICLAS_PATH/build_DSMC_2/bin/piclas

For the top boundary I use open boundary condition and

Part-Species1-Surfaceflux2-BC                   = 2            ! PartBound to be emitted from
Part-Species1-Surfaceflux2-velocityDistribution = maxwell_lpn
Part-Species1-Surfaceflux2-MWTemperatureIC      = 300
Part-Species1-Surfaceflux2-TempElec             = 300
Part-Species1-Surfaceflux2-Adaptive             = TRUE
Part-Species1-Surfaceflux2-Adaptive-Type        = 1
Part-Species1-Surfaceflux2-Adaptive-Pressure    = 0.2

There is still this issue:

Other entries from parameter.ini:

Part-Weight-SurfFluxSubSides = 200
AdaptiveBC-SamplingIteration                    = 100
AdaptiveBC-TruncateRunningAverage               = T       ! DEFAULT: F
Part-Weight-Type                = radial  
Part-Weight-Radial-ScaleFactor  = 600
Part-Weight-CloneMode           = 2
Part-Weight-CloneDelay          = 10

Do I need to do anything else in order to recompile piclas? The new build directory was empty, so I think the change in the file will be taken into account during compilation? Would this also work with existing build directories?

Let me know if I can provide any more details. Right now I am using the wall boundary as a workaround.

Best regards

Max

[pnizenkov]

Version 4.0.0 has been released and contains the mentioned bugfix. Let us know whether it fixes your issue!

[MaXXe4th]

Hey Paul, thanks for the information. I will try out as soon as possible and give you feedback.

Best regards

[MaXXe4th]

The first time steps were made with radial weighting and an adaptive boundary at the top and right edge of the domain:

! Radial Weighting
Part-Weight-Type                = radial    ! Activates a radial weighting in y for the axisymmetric simulation based on the particle position.
Part-Weight-Radial-ScaleFactor  = 20       ! Axisymmetric radial weighting factor, defining the linear increase of the weighting factor

Part-Species1-Surfaceflux3-BC                   = 6            ! PartBound to be emitted from
Part-Species1-Surfaceflux3-velocityDistribution = maxwell_lpn
Part-Species1-Surfaceflux3-MWTemperatureIC      = 300
Part-Species1-Surfaceflux2-PartDensity          = 4.83E+19
Part-Species1-Surfaceflux3-TempElec             = 300
Part-Species1-Surfaceflux3-VeloIC               = 0.0
Part-Species1-Surfaceflux3-VeloVecIC            = (/0.,1.,0./)
Part-Species1-Surfaceflux3-Adaptive             = TRUE
Part-Species1-Surfaceflux3-Adaptive-Type        = 1
Part-Species1-Surfaceflux3-Adaptive-Pressure    = 0.2

AdaptiveBC-SamplingIteration                    = 1000
AdaptiveBC-TruncateRunningAverage               = T       ! DEFAULT: F

t=0.005

Both boundaries seem to be working fine.

I then did a restart from the DSMC file in order to get sufficient particles per cell. (I aim for at least 10)

Particles-MacroscopicRestart = T
Particles-MacroscopicRestart-Filename        = DSMCResult.h5

[...]

Part-Weight-Type                    = cell_local
Part-Weight-CellLocal-QualityFactor = 0.8
Part-Weight-CellLocal-MinParticleNumber = 20
Part-Weight-CellLocal-MaxParticleNumber = 100
Part-Weight-CellLocal-SymAxis-MinPartNum = 50
Part-Weight-CellLocal-ApplyMedianFilter = T
Part-Weight-CellLocal-RefinementNumber = 2

The same boundary conditions now yield a reduced number density at the top boundary. The boundary at the right (separate boundary conditions in hopr.ini and parameter.ini) is set to exactly the same values.

t=0.0055

I attach the case files here, so you can check the settings: piclas_case.zip

Here are two more results for the restart from DSMCResult.h5:

sim_part_num for t=0.0055 weighting_factor for t=0.0055

[pnizenkov]

Hi Max, sorry for the delayed response! I could reproduce it with a smaller setup and it seems we have missed the same bug but for the case of an adaptive surface flux. You can try to add this code below line 1078 in particle_surface_flux.f90, this should fix your issue:

    ELSE
      ! No subdivision due to axisymmetry and weighting (for other particle weighting methods: weight is included in nVFR)
      CALL RANDOM_NUMBER(RandVal1)
      PartInsSubSide = INT(ElemPartDensity / Species(iSpec)%MacroParticleFactor * dtVar * nVFR + RandVal1)

I extended the test case to cover also the adaptive surface flux case and we will add the bugfix in the next release.

[MaXXe4th]

Thanks for the answer. I'll have a look soon and give you some feedback.

Best regards

Max