In a problem with a discontinuous internal boundary (a straight
converging-diverging nozzle) with a sharp throat, I get a NAN in the
reactive case when the detonation interacts with the throat. The
problem is however not encountered when the flow is non-reactive. Any
ideas about the culprit?
I have uploaded a script illustrating the problem called
DetonationCDRamp_EZ
matei
Sorry for being slow to reply, but I've been away on travel. Also, many happy returns on yesterday's anniversary!
On Wed, 26 Nov 2008, matei wrote:
> In a problem with a discontinuous internal boundary (a straight > converging-diverging nozzle) with a sharp throat, I get a NAN in the > reactive case when the detonation interacts with the throat. The > problem is however not encountered when the flow is non-reactive. Any > ideas about the culprit?
The inert case "working" says nothing about the reactive case. For you're using a linearized Riemann solver which is known to have problems in strong expansions, and so sooner or later you will generate NaN's.
Anyhow, I've just take a look at your script and you have a conceptual error to fix before you worry about the Riemann solver. The level set function used to define the surface is supposed to measure the distance from the surface. In your script you are using the DY from the surface and not the actual distance, which should be measured along the normal from a cell to the surface. Thus your inert case is actually incorrect.
I suggest you first fix the level-set function then we'll see what other issues rear their ugly head. You might also case to streamline you problem setup to take advantage of amr_sol's latest syntax: (i) use makefield {G$lmax}; (ii) replace def blocks by their fold::amr_sol counterparts.
I am still at a loss here with representing the level set
discontinuously. Is there any trick to
model the discontinuity in an acceptable way? although the normal is
not defined right at the corner?
Short of providing a level-set with a rounded corner, I have been
trying something on the lines of:
fold::amrita { define PSI
PSI1 ::= (((X[]-$X1)*($Y2-$Y1)-(Y[]-$Y1)*($X2-
$X1))) / sqrt(($X2-$X1)**2+($Y2-$Y1)**2)
PSI2 ::= -(((X[]-$X3)*($Y2-$Y3)-(Y[]-$Y3)*($X2-
$X3))) / sqrt(($X2-$X3)**2+($Y2-$Y3)**2)
wt12 ::= X[]-$X2 # cutoff to define the
level set piecewise
wt12 ::= wt12[]>0 ? 1 : 0
PSI ::= (1-wt12[])*PSI1[] + wt12[]*PSI2[]
#PSI ::= min(PSI1[],PSI2[])
}
AXS::InternalBoundary {
PSI = PSI[]
}
where the points 1, 2, 3 are the vertices of my CD nozzle.
The run starts out with a NAN, although the grid appears to be
correctly constructed. Any thoughts?
I have uploaded my script, it is called shockCD_new
matei
> Thanks for your input and pointing out my error.
> I am still at a loss here with representing the level set
> discontinuously. Is there any trick to
> model the discontinuity in an acceptable way? although the normal is
> not defined right at the corner?
You're thinking about things in the wrong way. The level set
gives the distance of a cell to the nearest solid boundary.
Thus in your example the level set in the vicinity of the throat
will be a series of arcs whose bounding sectors are the normals
to the front and back sides of the throat. Sketch this out by hand
and you'll see that it is a well-behaved function and there is
no discontinuity.
Moreover, since the normal used to impose the solid boundary function
is found by differentiating the level set, which is intrinsically smooth,
the normal is also well-behaved. The real weakness of the approach
has to do with how the flow data is reflected to achieve the desired
solid boundary. At a sharp-corner, such as the present throat, the "ghost" cell just inside the boundary receives data from two flow cells, the ones either side of the apex.
Anyhow, I've had a chance to investigate your script a bit further
and there is clearly something wrong with AXS::InternalBoundary.
For instance, try running the script with Y2 set to 7 and with lmax 0
and you'll see that a NAN is generated after 1 time step. It turns out that the implementation for the Aslam Xu Stewart
boundary consition, as given to me by Aslam, is weak in that
it falls over if a cell lies exactly on the boundary.
Your PSI1 and PSI2 are now correct.
The business with the wt12 you can forget about.
> AXS::InternalBoundary {
> PSI = PSI[]
> }
> where the points 1, 2, 3 are the vertices of my CD nozzle.
> The run starts out with a NAN, although the grid appears to be
> correctly constructed. Any thoughts?
With the above fixes you'll find that your script still blows up.
But now the failing is due to the fact that the linearized
Riemann solver is not positivity preserving and so pressure
can go negative. Consequently, unless you're
willing to implement a more robust solver, you have two options:
(i) increase the throat area, Y=7 works at this end; (ii)
try rounding the throat to reduce the strength of the expansion.
> I am still at a loss here with representing the level set
> discontinuously. Is there any trick to
> model the discontinuity in an acceptable way? although the normal is
> not defined right at the corner?
You're thinking about things in the wrong way. The level set
gives the distance of a cell to the nearest solid boundary.
Thus in your example the level set in the vicinity of the throat
will be a series of arcs whose bounding sectors are the normals
to the front and back sides of the throat. Sketch this out by hand
and you'll see that it is a well-behaved function and there is
no discontinuity.
[mir] Thanks. I'll admit I am a bit confused. What's the best
background reading to think of it the right way? Would the Xu, Aslam &
Stewart ref. be a good start to understand the implementation?
Moreover, since the normal used to impose the solid boundary function
is found by differentiating the level set, which is intrinsically
smooth,
the normal is also well-behaved. The real weakness of the approach
has to do with how the flow data is reflected to achieve the desired
solid boundary. At a sharp-corner, such as the present throat, the "ghost" cell just inside the boundary receives data from two flow cells, the ones either side of the apex.
Anyhow, I've had a chance to investigate your script a bit further
and there is clearly something wrong with AXS::InternalBoundary.
For instance, try running the script with Y2 set to 7 and with lmax 0
and you'll see that a NAN is generated after 1 time step. It turns out that the implementation for the Aslam Xu Stewart
boundary consition, as given to me by Aslam, is weak in that
it falls over if a cell lies exactly on the boundary.
Therefore you will need to modify the file:
$AMRITA/plugin/amr_sol/src/amrso/axs.F
Specifically, replace line 142:
DIST = PSI(I,J)
by:
DIST = PSI(I,J)
IF(DIST.EQ.0.0 D0) THEN
DIST = 1.0 D-10
ENDIF
then recompile amr_sol, using:
amrmake plugin::amr_sol
This fix is not particularly well thought out, but it
will do for now.
> Short of providing a level-set with a rounded corner, I have been
> trying something on the lines of:
Your PSI1 and PSI2 are now correct.
The business with the wt12 you can forget about.
> AXS::InternalBoundary {
> PSI = PSI[]
> }
> where the points 1, 2, 3 are the vertices of my CD nozzle.
> The run starts out with a NAN, although the grid appears to be
> correctly constructed. Any thoughts?
With the above fixes you'll find that your script still blows up.
But now the failing is due to the fact that the linearized
Riemann solver is not positivity preserving and so pressure
can go negative. Consequently, unless you're
willing to implement a more robust solver, you have two options:
(i) increase the throat area, Y=7 works at this end; (ii)
try rounding the throat to reduce the strength of the expansion.
> > I am still at a loss here with representing the level set > > discontinuously. Is there any trick to > > model the discontinuity in an acceptable way? although the normal is > > not defined right at the corner? > You're thinking about things in the wrong way. The level set > gives the distance of a cell to the nearest solid boundary. > Thus in your example the level set in the vicinity of the throat > will be a series of arcs whose bounding sectors are the normals > to the front and back sides of the throat. Sketch this out by hand > and you'll see that it is a well-behaved function and there is > no discontinuity. > [mir] Thanks. I'll admit I am a bit confused. What's the best > background reading to think of it the right way? Would the Xu, Aslam & > Stewart ref. be a good start to understand the implementation?
The Xu, Aslam, Stewart paper would be a good place to start. But there's no substitute for trying out lots of small examples. I would also suggest that you stick with the regular Euler equations when experimenting as then there are less things to go wrong. As always, you should look beyond the claims made in the formal write-up and view things with a critical eye.
> With the above fixes you'll find that your script still blows up.
> But now the failing is due to the fact that the linearized
> Riemann solver is not positivity preserving and so pressure
> can go negative. Consequently, unless you're
> willing to implement a more robust solver, you have two options:
> (i) increase the throat area, Y=7 works at this end; (ii)
> try rounding the throat to reduce the strength of the expansion.
There is another way, which I know James is not at all a fan of, but
which
is used in many schemes: modify the solver simply not to
let the pressure go negative, i.e. if P falls below P_tol, reset it to
P_tol,
and hope this only affects the solution locally. One would have to
check
the result is not sensitive to P_tol, provided is is small. It may or
may
not work in your case.
> > With the above fixes you'll find that your script still blows up. > > But now the failing is due to the fact that the linearized > > Riemann solver is not positivity preserving and so pressure > > can go negative. Consequently, unless you're > > willing to implement a more robust solver, you have two options: > > (i) increase the throat area, Y=7 works at this end; (ii) > > try rounding the throat to reduce the strength of the expansion.
> There is another way, which I know James is not at all a fan of, but > which > is used in many schemes: modify the solver simply not to > let the pressure go negative, i.e. if P falls below P_tol, reset it to > P_tol, > and hope this only affects the solution locally. One would have to > check > the result is not sensitive to P_tol, provided is is small. It may or > may > not work in your case.
---- ^ |
Much will depend on what your definition of "work" and "not sensitive" is.
But let me play devil's advocate for a second and ask: how many years do you expect to be performing reactive flow simulations? If the answer is more than 1 or 2 then you really ought to bite the bullet and become master of your own destiny. By that I mean you, or a student, should take an active interest in the innards of a patch-integrator. You'll be surprised that it is not much code at all and so there is really no excuse for not trying out algorithmic variations on a theme. And who knows, something might even filter back to the base AMRITA installation for others to benefit from.
In this case, Gary's suggestion has no real mileage in it. A better bet would be for you to read up on how Roe's scheme can be made positivity preserving by adjusting its wave speeds in a fashion that is rooted in numerical analysis rather numerical voodoo.
I did say James wouldn't like it...
However, I am not suggesting this as a permanent fix, just a third
option
in James' temporary fix list (i) and (ii) while the level set issues
are sorted out.
Of course, if one is going to fix this problem properly then the best
way to do it
is as James says have a solver "horse for the course" which can handle
it properly, and there would then be no need for James' (i),(ii) or my
(iii) suggestion.
Re: "numerical vodoo" - this is why something like Amrita is required
- many "modern" shock capturing schemes are full of all sorts of
fixes
and (potentially sensitive) numerical "tweak" parameters (often hidden
from the user) for "solving"
different problems. Additional artificial viscosity is another widely
used one to damp
numerical shock oscillations, etc, which is in widespread use, which
I'm would guess that
again James doesn't hold with(?). With Amrita, one can test the
effects on the solution of these
numerical fixes in a systematic way. Since, like it or lump, several
avialable schemes do
use the P_min "trick", with an Amrita-like facility one can actually
test quantitatively and
systematically what effect such fixes has on the solution for problems
with strong expansions, as compared to say positivity preserving
methods (and perhaps rule them out).
Gary
> > > willing to implement a more robust solver, you have two options:
> > > (i) increase the throat area, Y=7 works at this end; (ii)
> > There is another way, which I know James is not at all a fan of, but
> > which
> > is used in many schemes: modify the solver simply not to
> > let the pressure go negative, i.e. if P falls below P_tol, reset it to
> > P_tol,
> > and hope this only affects the solution locally. One would have to
> > check
> > the result is not sensitive to P_tol, provided is is small. It may or
> > may
> > not work in your case.
> ----
> ^
> |
> Much will depend on what your definition of "work" and "not sensitive" is.
> But let me play devil's advocate for a second and ask: how many years do
> you expect to be performing reactive flow simulations? If the answer is
> more than 1 or 2 then you really ought to bite the bullet and become
> master of your own destiny. By that I mean you, or a student, should take
> an active interest in the innards of a patch-integrator. You'll be
> surprised that it is not much code at all and so there is really no excuse
> for not trying out algorithmic variations on a theme. And who knows,
> something might even filter back to the base AMRITA installation
> for others to benefit from.
> In this case, Gary's suggestion has no real mileage in it. A better bet
> would be for you to read up on how Roe's scheme can be made positivity
> preserving by adjusting its wave speeds in a fashion that is rooted in
> numerical analysis rather numerical voodoo.
My likes and dislikes do not realy come into it. Matei could certainly bodge things today by fudging Pmin, but he would likely have to revisit the bodge for each new poblem he ran. And so all I was pointing out is that in the long-run it would likely be less effort to do things properly in the first place.
> However, I am not suggesting this as a permanent fix, just a third > option > in James' temporary fix list (i) and (ii) while the level set issues > are sorted out. > Of course, if one is going to fix this problem properly then the best > way to do it > is as James says have a solver "horse for the course" which can handle > it properly, and there would then be no need for James' (i),(ii) or my > (iii) suggestion. > Re: "numerical vodoo" - this is why something like Amrita is required > - many "modern" shock capturing schemes are full of all sorts of > fixes > and (potentially sensitive) numerical "tweak" parameters (often hidden > from the user) for "solving" > different problems. Additional artificial viscosity is another widely > used one to damp > numerical shock oscillations, etc, which is in widespread use, which > I'm would guess that > again James doesn't hold with(?). With Amrita, one can test the > effects on the solution of these > numerical fixes in a systematic way. Since, like it or lump, several > avialable schemes do > use the P_min "trick", with an Amrita-like facility one can actually > test quantitatively and > systematically what effect such fixes has on the solution for problems > with strong expansions, as compared to say positivity preserving > methods (and perhaps rule them out).
One could indeed do what you suggest. Does that mean you're volunteering to do so?
Anyhow, the relevant text from AMRITA's VKI notes is:
fold::latex { place Amrita in present context In the context of this lecture series, \Amrita\ can be viewed as a software system for automating CFD investigations; right down to the level of constructing documents which explain both the purpose and the outcome of a particular exercise. Automation is seem as the key to improving numerical reliability, repeatability and productivity to the point where algorithms could be improved through massed scrutiny. To reach this ideal, \Amrita\ strives to provide a computational framework which is sufficiently attractive to both novice and expert alike that it might help erode the present cottage-industry mentality, and its concomitant vagaries, where CFD codes are crafted on a one-off basis. Specifically, the many latitudes introduced by mundane activities such as preparing input files and post-processing results, together with the sheer drudgery of orchestrating investigations by hand, ensure that there is no convenient, watertight basis for the exchange of practical information to feed back and improve the underlying CFD algorithms. As a result, important observations can be slow to percolate through the research community.
}
the notes appeared February 1998, but the prose was actually written in 1996 as the forward to an "Introduction to Amrita" whose Chapter 2 first contained the my.script example.
After waiting 12+ years, it now seems that "reproducible research" is beginning to take off and so you might want to check out www.reproducibleresearch.org and www.reproducibility.org . The point I was making to Matei was that his long term plans impact on how he should be using a system like AMRITA. For example, consider the attached annotated listings, run_madagascar.pdf and ex5.pdf . The first is an AMRITA-ized version of a MADAGASCAR example from www.reproducibility.org and the second is an AMRITA-ized example of LeVeque's CLAWPACK.
It's a lot of effort to learn AMRITA's ropes if all you want to do is to run a black-box code to get out some numerical results. But given the scope of the system, the effort is actually quite small if you take a broader interest in its implications. For example, today someone might argue: why bother with AMRITA when we have Python, or Ruby? But such questions really miss AMRITA's point. As a case in point, the attached, co-operative-parsing.pdf shows that AMRITA is able to leverage off third-party software in nonstandard ways. But I suspect it will be another 12+ before the significance of re-entrant program folds is generally appreciated.
James
p.s. The attachements exceeded google's 4MB e-mail limit and so I'm sending them as separate messages.
Attached is an AMRITA-ized Clawpack example. It's just one example from a sequence, and you'll have to imagine it in the context of a conversation between two colleagues "Alice and Bob." If you've not come across this duo before, see:
> James Quirk wrote: > > Attached is an example, co-operative-parsing.pdf, that shows the > > nonsense of programming language advocacy.
> Neat.
> To be more idiomatic, the Ruby bits should be
> lines 11 & 33: puts "inside an enigma!"
> lines 42, 49, & 58: remove semicolons
> line 59: remove parentheses
The Ruby code was purposefully written in that form so that it would more closely match its Perl counterpart. But as the source lives inside the PDF, you could easly incorportate your suggested edits.
> > The Ruby code was purposefully written in that form > > so that it would more closely match its Perl counterpart.
> But Ruby != Perl ? It would seem to demonstrate a lack > of intimacy with the Ruby language to attempt to "closely > match its Perl counterpart".
I will gladly admit to a lack of intimacy with Ruby as I do not use it as often as I would like. And Ruby is clearly not Perl. But the, deliberate, thought I had in mind when I constructed the script was to show that there was direct overlap between programming languages. Obviously this runs the risk of being misunderstood. But when I can muster the energy, I will add a thread that addresses your concerns. Incidentally, there is a term for programs that are valid in more than one language, but it's slipped my mind.
> > But as the source lives inside the PDF, you could > > easly incorportate your suggested edits.
> Definitely, but as an example, it may be better received > by the various "camps" if each camp sees it's language > used in idiomatic fashion?
I do see where you're coming from. But quite frankly it's impossible to please all camps all the time, and I've long given up trying to do so.
In this instance the Perl-Ruby-Python example should have been taken in the context of the Clawpack and Madagascar examples that I also posted. You're worrying at the level of syntactic variations on a theme, which is your perogative. I'm worrying in a different direction. Namely, I believe the machinery needed to support reproducible research (RR) transcends niche technical fields. But given human nature, it's almost inevitable that in the coming years we'll see numerous RR efforts spring up, as different fields operate in splendid isolation from one another.
It would be nice to think that the many pronged approach would, via some form of genetic evolution, result in robust RR methods whose utility were widely appreciated. And this way well happen in the long term. But in the short term it will be disastrous, for there will be much duplication of effort and endless big-Endian/little-Endian arguments as how things should be done. For instance, once people twig to the power of a re-entrant fold-tree, there will inevitably be arguments as to what type of fold-markeres to use, and so on. Thus losing sight of the bigger scientific picture.
> They again, maybe I misunderstand the intent of the example.
Yes you have, but I thought I'd already put you straight when we spoke on the phone?
|
