Determining the cause of a numerical instability

[Bjoern Bingham]

I am running a case where the only difference from jobs that ran with no errors is the gravitational constant.  I am simulating smoke transport under Earth, Martian, Lunar, and Microgravity. This was the Micro (zero) gravity case. I am running a locally compiled FDS on a cluster:

Fire Dynamics Simulator

 Current Date     : September 29, 2023  09:02:12
 Revision         : FDS-6.8.0-406-gc36eed0-master
 Revision Date    : Thu Jul 20 04:43:12 2023 -0400
 Compiler         : ifort version 2021.8.0
 Compilation Date : Jul 21, 2023 12:22:55

 MPI Enabled;    Number of MPI Processes:      27
 OpenMP Enabled; Number of OpenMP Threads:      2

 MPI version: 3.1
 MPI library version: Intel(R) MPI Library 2021.8 for Linux* OS

Approximately 45 seconds of the intended 600, into the simulation, it fails and the only information is a Numerical Instability with no particular pointers as to why. The smokeview file will not open so I cannot check for a root cause that way. Can anyone instruct me on where to look for information on why it would have failed?  I supply the input file or output files if that would help

I really appreciate any help you can provide.

Bjoern B

[Randy McDermott]

Sending the input file would help.  Thanks

--
You received this message because you are subscribed to the Google Groups "FDS and Smokeview Discussions" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fds-smv+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/fds-smv/e755d460-248d-464b-aba7-f3993cfb4e3fn%40googlegroups.com.

[Bjoern Bingham]

The input file is attached.  Had to zip it to attach. 

[Kevin McGrattan]

I am running the case.

[Kevin McGrattan]

The instability occurs just above the little purple  vent, but I do not know exactly why. I ran the case with just the lower, central mesh and it ran fine.  I'll keep looking. 

One suggestion is to simplify the case more and more. It is hard to detect the problem with so many lines of input. 

Has this case ever run successfully? If so, what is different now?

[Jonathan Hodges]

You might want to make the exterior mesh boundaries OPEN as well, or at least one of them. The "interior wall" surface has an exposed backside boundary condition so you may be adding heat to the exterior. Since that region is fully enclosed in pressure zone 1 without any release you could get a numerical instability over time building up there. You also have inlet and exhaust fans which are both specified as volume flow rates. If the temperature of the extracted gas is elevated the mass flow will be different between your inlets and outlets. Since you have an open vent that should not lead to a numerical instability, but it is something worth mentioning.

[Kevin McGrattan]

I had that same thought. I'm running the case now.

[Kevin McGrattan]

Spoke too soon. It failed again.

[Kevin McGrattan]

My next suspect are these 1 micron particles. 

[Bjoern Bingham]

Yeah, that was an issue prior to adding the open vent. The only difference between this case and others that have all run successfully is GVEC = 0.0 rather than a nonzero value. I can try to simplify the case. I added quite a few devices and switched to a heat flux rather than a set ramped surface temperature at the small purple vent. I also defined a size distribution for the particulate. Prior to these changes, I had microgravity cases run successfully.

[Jonathan Hodges]

Here is a version where I removed all the OBST and replaced with a convex hull of the original shell.

[Kevin McGrattan]

I removed the particles and the case is running past its prior point of failure.  Tiny particles can be tricky because FDS treats them as independent individual particles (i.e. water droplets) rather than a cloud of dust.  Assuming the case runs without particles, you might try reducing the drag coefficient. With such tiny particles, it's hard to say what the drag coefficient is. The particles just go with the flow. You might also try launching more particles per second to avoid having too much mass concentrated on a single particle. 

[Kevin McGrattan]

I have tested a few ideas on this case. My latest test involves the heat transfer coefficient to the tiny sphere. We use a correlation for the Nusselt number of a sphere. Then the HTC is

h = Nu k / L

L is the length scale, which in your case is the diameter of this tiny particle. That makes h huge, and this might lead to a numerical instability because the temperature of the particle changes too much in a single time step. Two options I am testing:

1. CHECK_HT=T will reduce the time step to account for the ;large HTC. 

2. HEAT_TRANSFER_COEFFICIENT=50 (a guess) on the SURF line, which will limit the heat transfer to the particle in a single time step.

In reality, these particles are so small that they take on the temperature of the gas. I'll let you know.

[Kevin McGrattan]

Setting HTC=50 did not prevent the instability. I am now running the case with the following changes from the original file:

< &SURF ID = 'PTFES', MATL_ID = 'PTFE', THICKNESS = .00000001, GEOMETRY='SPHERICAL'/
---
> &SURF ID = 'PTFES', MATL_ID = 'PTFE', THICKNESS = .00000001, GEOMETRY='SPHERICAL', HEAT_TRANSFER_COEFFICIENT=50 /
23c23
< &PART ID = 'PMP', MAXIMUM_DIAMETER = 1.0, MINIMUM_DIAMETER = 0.01, SURF_ID = 'PTFES'  , CHECK_DISTRIBUTION = T, QUANTITIES = 'PARTICLE AGE', 'PARTICLE DIAMETER', 'PARTICLE MASS', MONODISPERSE = F,  CNF_RAMP_ID = 'dist', N_STRATA = 10/
---
> &PART ID = 'PMP', MAXIMUM_DIAMETER = 1.0, MINIMUM_DIAMETER = 0.01, SURF_ID = 'PTFES'  , CHECK_DISTRIBUTION = T, QUANTITIES = 'PARTICLE AGE', 'PARTICLE DIAMETER', 'PARTICLE MASS', 'PARTICLE TEMPERATURE', MONODISPERSE = F,  CNF_RAMP_ID = 'dist', N_STRATA = 10, SAMPLING_FACTOR=1, DRAG_COEFFICIENT=0.4 /
309c309
<       NET_HEAT_FLUX=15.224 /
---
>       NET_HEAT_FLUX=15.224,NPPC=10, DT_INSERT=0.005 /

HTC=50 prevents the particles from heating up too fast over a single time step.

DRAG_COEFFICIENT=0.4 prevents the particle drag force from being a huge number because of the tiny diameter and low Reynolds number. 

NPPC=10 tells FDS to launch more particles from the surface. This may not affect the case significantly, but it will (along with SAMPLING_FACTOR=1) show more particles in Smokeview.

I have also launched a case with CHECK_HT=T, but that case is still running. It has the effect of significantly lowering the time step, which sometimes overcomes these kinds of instabilities. 

[Kevin McGrattan]

CHECK_HT lowered the time step by a factor of 10, but the case still failed. For the moment, I think setting the DRAG_COEFFICIENT to a reasonably small number will work for you. Essentially, the particles go with the flow.

[Bjoern Bingham]

Thanks for continuing to look into this Kevin. I will try setting the DRAG_COEFFICIENT to 1 to start with. Based on some quick math, assuming a small difference in particle and gas velocity and a mode particle size of 100nm, my DC should be closer to 2000. I will keep this thread updated with my progress.

[Bjoern Bingham]

With setting DRAG_COEFFICIENT =1, the simulation made it past its original failure point but still crashed ~270 seconds in, I will lower the DC more and try again.

[Kevin McGrattan]

For me it ran to completion with DC=0.4. That is the lower bound of the drag coefficient.

[Bjoern Bingham]

Ok, thank you, I will try that.

[Bjoern Bingham]

Hi Kevin, 

Unfortunately, with DC = 0.4, the case is failing closer to its original failure point. I have tried a couple of other drag coefficients between 1 and 0.4 but with no success. I will continue trying slightly higher DCs to see if that improves the simulation performance.

Apologies for the delayed update. I spent the remainder of last week focusing on other tasks.

[Bjoern Bingham]

Is DC = 0.4 the only change you made to my original file?

[Kevin McGrattan]

I will check. I am not in the office now.

--

You received this message because you are subscribed to the Google Groups "FDS and Smokeview Discussions" group.
To unsubscribe from this group and stop receiving emails from it, send an email to fds-smv+u...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/fds-smv/18fde63c-1ac7-42f6-b678-88babd22c730n%40googlegroups.com.

[Kevin McGrattan]

[mcgratta@burn Test]$ diff fine_MicroG_HFlux_mpi_929.fds fine3.fds
8c8
< &HEAD CHID='fine_MicroG_HFlux_mpi_929' TITLE='This is an example case' /
---
> &HEAD CHID='fine3' TITLE='This is an example case' /
21c21

< &SURF ID = 'PTFES', MATL_ID = 'PTFE', THICKNESS = .00000001, GEOMETRY='SPHERICAL'/
---
> &SURF ID = 'PTFES', MATL_ID = 'PTFE', THICKNESS = .00000001, GEOMETRY='SPHERICAL', HEAT_TRANSFER_COEFFICIENT=50 /
23c23
< &PART ID = 'PMP', MAXIMUM_DIAMETER = 1.0, MINIMUM_DIAMETER = 0.01, SURF_ID = 'PTFES'  , CHECK_DISTRIBUTION = T, QUANTITIES = 'PARTICLE AGE', 'PARTICLE DIAMETER', 'PARTICLE MASS', MONODISPERSE = F,  CNF_RAMP_ID = 'dist', N_STRATA = 10/
---
> &PART ID = 'PMP', MAXIMUM_DIAMETER = 1.0, MINIMUM_DIAMETER = 0.01, SURF_ID = 'PTFES'  , CHECK_DISTRIBUTION = T, QUANTITIES = 'PARTICLE AGE', 'PARTICLE DIAMETER', 'PARTICLE MASS', 'PARTICLE TEMPERATURE', MONODISPERSE = F,  CNF_RAMP_ID = 'dist', N_STRATA = 10, SAMPLING_FACTOR=1, DRAG_COEFFICIENT=0.4 /
309c309
<       NET_HEAT_FLUX=15.224 /
---
>       NET_HEAT_FLUX=15.224,NPPC=10, DT_INSERT=0.005 /

 Fire Dynamics Simulator

 Current Date     : October  2, 2023  10:27:04
 Revision         : FDS-6.8.0-635-g2e1bf02-master
 Revision Date    : Mon Sep 18 10:20:24 2023 -0400
 Compiler         : ifort version 2021.7.1
 Compilation Date : Sep 18, 2023 15:48:02

 MPI Enabled;    Number of MPI Processes:      27

 OpenMP Disabled

 MPI version: 3.1
 MPI library version: Intel(R) MPI Library 2021.7 for Linux* OS

[Bjoern Bingham]

Ok, thank you, Kevin.

I will make the changes and hopefully, that fixes my issue.

[Bjoern Bingham]

This worked.

Is there any reason why changing the vent locations would cause this not to work? In the original case, which has now run successfully with the changes, the ventilation was top-down. I am also comparing flow orientations and the across-flow case is still failing due to an unidentified numerical instability.

[Kevin McGrattan]

I would think it is still an issue with the tiny particles. Can you post the case that fails.

[Bjoern Bingham]

Here it is

[Kevin McGrattan]

I decreased the DRAG_COEFFICIENT to 0.04 and the case has passed its point of failure. I think the problem is related to assumptions that we are making about the drag from individual particles summed within a single grid cell. Ideally, as the particles get smaller and smaller, they should just go with the flow. Somehow that is not happening here, so I am lowering the drag coefficient to avoid some overshoot in the transfer of momentum from particles to gas. I'll keep watching it. 

[Kevin McGrattan]

The case with the reduced drag coefficient is still running. I think this solution will work for you temporarily until we improve how we handle these tiny particles.

[Bjoern Bingham]

Ok, thank you, Kevin. I appreciate your effort on this.