Opened 10 years ago

Closed 9 years ago

Last modified 5 years ago

#67 closed defect - convergence (fixed)

failed assert in rt_diffuse.cpp at line number 393

Reported by: peter Owned by: peter
Priority: major Milestone: c13 release
Component: thermal convergence Version: trunk
Keywords: Cc:

Description (last modified by peter)

The attached test case crashes the mainline when modeling a high excitation PN with a 1/r2 density law outside the ionized region. It uses the DLAW command and the dense_fabden.cpp file is attached as well. An analysis by Ryan Porter yielded the following:

Here's where the problem stands then. At line 393 of rt_diffuse the below assert values. The comment accompanying it may well provide the proper hint:

/* nflux was reset upward in ConvInitSolution to encompass possible line and continuum emission. this test should not possibly fail. It could if the ionization were to increase with depth This is important because the nflux cell in ConInterOut is used to carry out the unit integration, and if it gets clobbered by diffuse emission the code will declare insanity in PrtComment */

ASSERT( ip >= 0 && ip < rfield.nflux );

When the assert is thrown ip == 3146, rfield.nflux == 3141, nelem == 25, and ion == 16.

Attachments (2) (1022 bytes) - added by peter 10 years ago.
input script
dense_fabden.cpp (511 bytes) - added by peter 10 years ago.
density law

Download all attachments as: .zip

Change History (11)

Changed 10 years ago by peter

Attachment: added

input script

Changed 10 years ago by peter

Attachment: dense_fabden.cpp added

density law

comment:1 Changed 10 years ago by peter

Description: modified (diff)

Problem was reported by Patrizia Manzitto.

comment:2 Changed 9 years ago by peter

Owner: changed from nobody to peter
Status: newassigned

On the trunk at r2716 this sim now crashes on a series of messages like this:

 PROBLEM DISASTER - the kinetic temperature, 2.771e+00K, is below the lower limit of the code, 2.800e+00K.
 This calculation is aborting.
 PROBLEM DISASTER - the kinetic temperature, 2.716e+00K, is below the lower limit of the code, 2.800e+00K.
 This calculation is aborting.
    ....... etc ........

This is caused by a failure of the temperature solver to bracket the zero point of the heating-cooling curve, causing it to assume that it reached a thermal front and taking a deep plunge. However, it never finds a solution at lower temperatures, so it eventually drops below 2.8 K and the code aborts.

comment:3 Changed 9 years ago by peter

Further analysis shows that the code correctly detects that the heating and cooling curves detach and that it needs to bifurcate to find a new solution. The current logic briefly searches at higher temperatures (this is needed for regions where the temperature is slowly rising) and then starts a long hunt downwards in the assumption that the new equilibrium must be at much lower temperatures. Normally that is the correct choice, but not in this case. Here the new equilibrium is at higher temperatures. To be more precise, the code should have taken a jump up from 17150 K to 35050 K. This is most likely due to the 1/r2 drop-off of the density. The temperature would then continually rise to 96590 K and then another bifurcation would occur to 2.576 MK.

Note that the stopping criterion in this sim is incorrect. It is set to stop at 100 K, but the temperature never drops to values near that temperature and keeps rising going out. At the point where the first bifurcation occurs, the density has already dropped to 2.113e-05. That is well below ISM values. If we would fix this particular problem, the sim would continue integrating outwards with ever rising temperatures until eventually it exceeds 10 GK and the code crashes again. We cannot fix that problem.

comment:4 Changed 9 years ago by peter

With c08_branch at r2695, the code integrates outwards until it reaches temperatures of 55.58 MK and then aborts with the following message:

PROBLEM DISASTER the density of ion 3 of element Lithium    is too small to be computed on this cpu.
Turn off the element with the command ELEMENT XXX OFF or do not consider gas with low density, the current hydrogen density is 4.55e-08 cm-3.

This is well below the lowest density we treat, so I would say that this is a wontfix for c08_branch.

comment:5 Changed 9 years ago by peter

Component: radiative transferthermal convergence
Type: defect - failed assertdefect - convergence

The original problem was fixed in r2045 on the mainline and r2047 on c08_branch.

comment:6 Changed 9 years ago by peter

Milestone: C10 branch
Version: trunk

PR #102 is a duplicate of this PR.

comment:7 Changed 9 years ago by peter

Resolution: fixed
Status: assignedclosed

The second problem is fixed on the mainline in r3349.

comment:8 Changed 9 years ago by peter

The test case slow/ also has a big thermal front to higher temperatures in the last zone. It was originally included in PR #102, but has been removed there since it was not related to the other problem that was reported there. The failure in slow/ was also solved by r3349.

comment:9 Changed 5 years ago by Ryan Porter

Milestone: C10 branchc13 release

Milestone C10 branch deleted

Note: See TracTickets for help on using tickets.