problem with discontinuous level set boundary in the reactive case
matei
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
James Quirk
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.
James
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?
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
James Quirk
On Wed, 3 Dec 2008, matei wrote:
>
> James,
>
> 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?
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.
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:
>
> 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[])
> }
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?
>
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.
James
Matei Radulescu
>
> 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.
background reading to think of it the right way? Would the Xu, Aslam &
Stewart ref. be a good start to understand the implementation?
James Quirk
>
>
> >
> > 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.
James
Gary Sharpe
> 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.
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.
James Quirk
----
^
|
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.
James
Gary Sharpe
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)
James Quirk
On Sat, 6 Dec 2008, Gary Sharpe wrote:
>
> I did say James wouldn't like it...
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.
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.
James Quirk
James Quirk
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:
http://en.wikipedia.org/wiki/Alice_and_Bob
James
Bil Kleb
> 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
Regards,
--
Bil Kleb
James Quirk
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.
James
>
> Regards,
> --
> Bil Kleb
>
Bil Kleb
> Bil,
Hi.
> 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".
> 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?
They again, maybe I misunderstand the intent of the example.
Regards,
--
Bil Kleb
http://fun3d.larc.nasa.gov
James Quirk
On Tue, 16 Dec 2008, Bil Kleb wrote:
>
> James Quirk wrote:
> > Bil,
>
> Hi.
>
> > 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?
James