Numerical instability

[gforney]

Application Version: 5.4.3 Serial
SVN Revision Number: 5540
Compile Date: 5 february 2010
Operating System: Windows

Compiled svn 5540 with GNU compilers and performed a calculation of the
attached file. After 880.1 seconds the calcualtaion is terminated due to
numerical error.

Tried to use smokeview to see if it could be possible to "see" the error
but I could not see anything. Do not know what is wrong so I have not
slimmed down the input file so that the error occurs quickly.

Since I am interested in heat transfer to surfaces I tried to set H_EDDY =
.TRUE. (or only T) but then the numerical instability occurs after 1.69
seconds instead. (the same result with svn 5587)

Joakim
(I know the file is very similar to the SP fire test in the validation suite)

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 12:02:49

---

• _Attachment: Rome_38.fds_

[gforney]

Added H_EDDY = T and CSMAG in the MISC line. These where the results (so far)

CSMAG = 0.17   => Numerical instability after 1.27 s
CSMAG = 0.20   => Numerical instability after 1.69 s
CSMAG = 0.25   => Still running after 16.77 s

Do not know if they help.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 13:36:00

[gforney]

H_EDDY is a new feature that we are still testing. What happens when you set it 
to .FALSE.? Also, increasing or decreasing CSMAG just increases or decreases the 
amount of numerical viscosity. By increasing to 0.25, you are effectively damping 
out the instability, but you are also inhibiting eddy formation. Let's see if H_EDDY

is the problem.

Original issue reported on code.google.com by mcgratta on 2010-02-15 13:53:44

[gforney]

I assume that H_EDDY is set to .FALSE. by default. In that case the numerical
instability will occur after 880.1 seconds. I'll run a calculation just to confirm,
it just takes a while...

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 15:10:00

[gforney]

Just checked the last of the three above.

CSMAG = 0.25   => Numerical instability after 23.96 s

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 15:11:34

[gforney]

Try setting BAROCLINIC2=F on MISC.  We are also testing a new baroclinic correction
and this may be in play.  Do not play with CSMAG.  Leave it as the default.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 15:16:45

• Labels added: Milestone-FDS_5.5

[gforney]

Note the "2" at the end of BAROCLINIC2!!

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 15:17:32

[gforney]

Would this be a correct MISC-line to try?

&MISC SURF_DEFAULT       = 'LWC WALL'
      TMPA               = 20.
      H_EDDY             = .TRUE.
      BAROCLINIC2        = .FALSE.
      RESTART            = .FALSE. /

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 15:31:48

[gforney]

Yes.  But let's also decrease the grid resolution and see if we can get this thing to
fail faster.  I am working with half the resolution right now.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 15:34:27

[gforney]

As you said, it will take a while to see if we get past the 880 sec barrier.

If it does not, you may want to try FDS6.  I recommend trying the following at the
lower resolution:

&MISC FDS6=T / (no need to touch H_EDDY or BAROCLINIC2)
&REAC HRRPUA_SHEET=1.E10, C_EDC=1./

In pool fire cases we have found that the HRRPUA_SHEET bound is controlling the
physics and it is better to get rid of it (this bound was implemented before the Eddy
Dissipation Concept model was in use).  But that leaves C_EDC unvalidated with the
new set of parameters.  So, this is the knob to turn to see if you can get your case
to match the data.  C_EDC=1 is a good place to start.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 15:49:05

[gforney]

With half the resolution and this MISC-line

&MISC SURF_DEFAULT       = 'LWC WALL'
      TMPA               = 20.
      H_EDDY             = .TRUE.
      BAROCLINIC2        = .FALSE.
      RESTART            = .FALSE. /

gave numerical instability after 2.29 s. 

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 15:49:14

[gforney]

I'll try the above and let it run tonight, my workday is almost over now and I cannot
check the calcualtions from home.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 15:51:56

[gforney]

Okay.  Let me try that scenario again.  I thought that one had worked for me.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 15:52:04

[gforney]

Ugh.  Here's where things get hard, because that case works for me.  Here are my mesh
lines...

&MESH   IJK = 24,24,24
    XB  = 0.0, 2.4, 0.0, 2.25, 0.0, 2.4 /

&MESH   IJK = 24,24,24
    XB  = 0.0, 2.4, 2.25, 4.5, 0.0, 2.4 /

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 15:55:21

[gforney]

&MISC SURF_DEFAULT       = 'LWC WALL'
      TMPA               = 20.
      H_EDDY             = .TRUE.
      BAROCLINIC2        = .TRUE.
      RESTART            = .FALSE. /

Numerical instability after 11.05 s

Get back with the rest tomorrow.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-15 15:55:28

[gforney]

Joakim,

To take the compiler issues off the table, here is the executable I am using:

http://fds-smv.googlecode.com/files/fds_dev.exe

R

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 16:04:36

[gforney]

I got a crash at 111 sec with FDS6=T.  I'm going to try a single mesh case.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-15 17:29:29

[gforney]

&MISC SURF_DEFAULT       = 'LWC WALL'
      TMPA               = 20.
      FDS6               = .TRUE.
      RESTART            = .FALSE. /

&REAC HRRPUA_SHEET = 1.E10
      C_EDC        = 1. /

Double the grid size   => Crashed at 130.49 sec.
"Normal" grid size     => Crashed at 94.47 sec.

I'll do the same runs today but with your executable.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-16 07:16:30

[gforney]

For information; my mesh for the "double the grid size" is:

&MESH   IJK = 24, 45, 24
    XB  = 0.0, 2.4, 0.0, 4.5, 0.0, 2.4 /

I now have all the previous crashes running with your executable.

The OS are one XP and the other Windows 7. The first crashes occurred on the XP
machine but a few of the latest have been on the W7 machine. (It is a lot faster...)

Original issue reported on code.google.com by jsandstrom42 on 2010-02-16 07:37:39

[gforney]

Double the grid size crashed at 118.28 sec with fds_dev.

I slimmed down the input file to get rid of all the extra output that is not needed
right now.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-16 09:07:29

---

• _Attachment: Rome_slim.fds_

[gforney]

Okay.  Thanks for your efforts.  I'll work on this today.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-16 12:57:48

[gforney]

Just let me know if there is something I can do to help.

&MISC SURF_DEFAULT       = 'LWC WALL'
      TMPA               = 20.
      H_EDDY             = .FALSE.
      BAROCLINIC2        = .TRUE.
      RESTART            = .FALSE. /

is still running after 1750 sec. (coarse grid)

Original issue reported on code.google.com by jsandstrom42 on 2010-02-16 13:41:22

[gforney]

The above calcuation ran through to the end.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-16 15:09:32

[gforney]

Joakim,

In FDS5, I think the main problem was being caused by the wall damping factor.  I
removed that for SVN 5601.  BAROCLINIC2=T is default for this version, so you will
only need this MISC line.

&MISC SURF_DEFAULT       = 'LWC WALL'
      TMPA               = 20.
      RESTART            = .FALSE.
      H_EDDY             = .TRUE. /

Unfortunately, this creates a problem for the H_EDDY model.  Quite simply, we need
some verification tests to see if this model is working correctly.  But since we are
focusing on FDS6 development it does not make much sense to put in the machinery
necessary for the damping factor.  I would note that Moinuddin and Li's
implementation is not complete because they only add a damping factor (not correctly
computed) for the heat transfer coefficient, but not for the viscosity near the wall
in general.  This issue is handled more gracefully by the dynamic model.  Further,
for FDS6 we will likely go with a log law model for the convective boundary condition
as is employed by atmospheric models.  So, again, it does not make much sense to
modify the current code to develop H_EDDY further.

I did still get an instability with FDS6, but I think this has more to do with the
REAC parameters we are testing.  Now that I think of it, I will take H_EDDY off of
FDS6 default until this is resolved.

So, let me know if things now work for you.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-17 00:00:34

[gforney]

Is it possible to test the log law model within the FDS5 or do I have to test it with
the rest of the FDS6 code? 

I am very interested in testing different convective heat transfer models and their
impact on the adiabatic surface temperature. There can be a difference between
different models since differences in h can have a large impact on the resulting AST.
Especially in early phases of a fire but also in the later stage.

Tested SVN 5602 with the above input for MISC and got numerical instability after
49.24 sec.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-17 11:07:11

[gforney]

I have not implemented log law model yet.  I'll try to get to that today.  No time
like the present...

But it is still curious that H_EDDY failed for you. This case worked for me in FDS5
but not FDS6.  But again, I think that has to do with REAC parameters, not h.  

Original issue reported on code.google.com by randy.mcdermott on 2010-02-17 13:07:52

[gforney]

Oh, try removing the REAC line you added so that you are just using FDS5 defaults,
except H_EDDY.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-17 13:09:33

[gforney]

I'll try that. Get back with the results tomorrow.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-17 16:23:14

[gforney]

Joakim,

I think I actually made some meaningful progress on this.  Several of the things I
was playing with before turned out to be red herrings.  I think the main issue was
the stability check based on the divergence.  In short, I modified the code (SVN
5604) to include the divergence TOGETHER with the CFL, instead of separately.  There
is a good theoretical reason behind this, which I'll explain in a moment.

But first note that I added the crude damping factor back to H_EDDY.  I also started
to code up the log law model, but it is even more clear now that we need to bite the
bullet and store the wall friction velocity in a wall cell array.  This is needed for
H_EDDY, the log law model, and Jason can use it for soot deposition.  So, stay tuned
and as I get that implemented I will bring the log law model on line too.  In case
you are interested, the paper I am working from for the log law is: Stoll, R.,
Porte-Agel, F. (2008) Large-Eddy Simulation of the Stable Atmospheric Boundary Layer
using Dynamic Models with Different Averaging Schemes. Boundary-Layer Meteorology,
126:1-28.

I have attached the input file I currently have running (+238 sec).  Note that I put
back the REAC line.  I want to encourage Kevin to make this change for the next
Validation Suite run (Kevin, we should talk off line).

So, let me know how this works for you.  I appreciate your help in working on this
Issue!

Randy

p.s. The theory behind the div constraint is as follows.  An explicit Euler update
of
the density equation is:

rho^{n+1} = rho^n - ( u^n * grad(rho)^n + rho^n * div(u)^n )

Therefore, the constraint that rho^{n+1} > 0 leads to the following time step
restriction:

dt < rho^n/( u^n * grad(rho)^n + rho^n * div(u)^n )

We can argue that the case we are most concerned with is when rho^n is small.  A
reasonable approximation to the denominator then is (superscript dropped):

dt < rho/( u*(rho-0)/dx + rho*div(u) )

or

dt < 1/( u/dx + div(u) )

The first part is the standard CFL constraint.  But this analysis says we should
consider the divergence together with the CFL.  Before we had been treating them
individually, and this made more of a difference than we expected.

Original issue reported on code.google.com by randy.mcdermott on 2010-02-17 19:41:58

---

• _Attachment: Rome_slim.fds_

[gforney]

Joakim,

Still having problems with FDS6.  However, I made some progress on H_EDDY.  See Issue
971.

R

Original issue reported on code.google.com by randy.mcdermott on 2010-02-18 00:39:06

[gforney]

I am happy to help. This would be my way to contribute to the project and your great
work.

I started a calculation with H_EDDY and the REAC as above with SVN 5608 (the latest).
I took the original file I submitted and am very excited to see if it works. It will
take a while though but I hope it will be worth it. I'll get back.

I'll happily try out the log-law model for comparison when it's implemented. Thanks
for the paper, I appreciate it.

Original issue reported on code.google.com by jsandstrom42 on 2010-02-18 08:55:10

[gforney]

I am closing this issue.  If there is an issue with latest verification or validation
case we can reopen.

Original issue reported on code.google.com by randy.mcdermott on 2011-08-26 19:40:13