<div dir="ltr">Hi Norbert,<div><br></div><div>I agree that the internal potential formula, Eqn. 25.39, reduces to the proper value if A1 and A3 both go to 2/3.  The issue is that in Simulation_init.F90 AA3 is used to calculate pressure:</div>
<div><br></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div><div>  p0     = (2.0/3.0)*sim_pi*sim_Newton*sim_density**2*amean**2</div></div><div><div>  sim_Pconst = p0 * AA3 * (1.0-sim_eccentricity**2)**(2.0/3.0)</div>
</div></blockquote><div><br></div><div>If we set eccentricity to zero, this reduces to the proper sim_Pconst = 2/3*sim_pi*sim_Newton*sim_density**2*amean**2 central pressure for a spherical object of constant density, but that's only because AA3 is artificially set to 1 when sim_eccentricity = 0.  At sim_eccentricity << 1, but > 0, sim_Pconst is incorrect since AA3 goes to 2/3 rather than 1.</div>
<div><br></div><div>The reason I'm running these tests is because I'm running post-merger evolution for merged binary white dwarf systems, and I'm seeing substantial differences between FLASH, SPH and moving mesh codes.  The remnant has enough similarities to either a MacLaurin or Jacobi spheroid that I'm hoping that I can spot the same differences by running spheroids through all the codes.  Since I've set gamma = 10, I don't think the initial pressure distribution matters a whole lot - the star doesn't redistribute mass very severely, so this bug, if it's a bug, shouldn't be an enormous deal.</div>
<div><br></div><div>I'm a bit concerned about approximating an incompressible fluid though.  Bob Fisher's already warned me that v/c_s << 1 systems tend to take many more timesteps and accumulate integration errors as a result.  Are there other particularly nasty known issues related with setting the EOS gamma to an unusually large value?</div>
<div><br></div><div>Thanks again for your help,</div><div><br></div><div>Charles</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Aug 15, 2013 at 11:07 AM, Norbert Flocke <span dir="ltr"><<a href="mailto:flocke@flash.uchicago.edu" target="_blank">flocke@flash.uchicago.edu</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Charles,<br>
<br>
You are right, I did not consider enough terms. On the other hand, the<br>
AA3 going to 2/3 with e -> 0 is the correct behavior. Please look at the<br>
FLASH user's guide in the Gravity Test Problems / MacLaurin (Section 25.3.4). There you will see the formula for the potential in terms of<br>
two constants A1 and A3. When e -> 0 both go to 2/3 and the result is the potential of a uniform density sphere.<br>
<br>
The MacLaurin main simulation was written a long time ago and is somewhat hard to read. I reimplemented the MacLaurin test problem as a unit test for the multipole gravity section (this sits in /SimulationMain/unitTest/<u></u>Multipole). Please have a look at the routine:<br>

<br>
   sim_initBlockAnalytical.F90<br>
<br>
where all these constants are properly defined and their correct behavior is listed. If you still want to use the old MacLaurin code, then you would have to go through that code and see what it exactly does. What seems strange to me at first sight is that AA3 is set to 1.0 when e -> 0. But maybe somewhere down the path this is the correct setting.<div class="HOEnZb">
<div class="h5"><br>
<br>
Cheers,<br>
Norbert<br>
<br>
<br>
<br>
<br>
On Wed, 14 Aug 2013, Chenchong Zhu wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Norbert,<br>
<br>
Thanks for the response!<br>
<br>
For sqrt(1-e^2)*asin (e), I don't think you kept enough of the expansion:<br>
<br>
asin(e) -> e + e^3/6<br>
sqrt(1 - e^2) -> 1 - e^2/2 - e^4/8 -> 1 - e^2/2<br>
<br>
AA3 = 2/e^2 - 2*sqrt(1-e^2)*asin (e) / e^3<br>
      = 2/e^2 - (2/e^3)*(e - e^3/3 - e^5/12)<br>
      = 2/e^2 - 2/e^2 + 2/3<br>
      = 2/3<br>
<br>
Throwing the expression into Python and setting e = tiny number gives the<br>
same result.<br>
<br>
Best,<br>
<br>
Charles<br>
<br>
<br>
On Wed, Aug 14, 2013 at 1:29 PM, Norbert Flocke<br>
<<a href="mailto:flocke@flash.uchicago.edu" target="_blank">flocke@flash.uchicago.edu</a>><u></u>wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Sorry, there is a typo in the last AA3 equations. It should read:<br>
<br>
<br>
  AA3 -> 2/e^2 - 2*(1-e^2/2)*e / e^3<br>
      -> 2/e^2 - 2/e^2 + (2*e^2/2*e^2)<br>
      -> 1<br>
<br>
Norbert<br>
<br>
<br>
On Wed, 14 Aug 2013, Norbert Flocke wrote:<br>
<br>
 Hi Charles,<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
The answer to your first question I can provide. As written in the code<br>
on line 78, the value of AA3 is<br>
<br>
 AA3 = (2*sqrt(1-e^2)/e^2) * (1/sqrt(1-e^2) - asin (e)/e)<br>
<br>
or, written out<br>
<br>
 AA3 = 2/e^2 - 2*sqrt(1-e^2)*asin (e) / e^3<br>
<br>
where e is the eccentricity. If e -> zero, then the MacLaurin (what<br>
coincidence!) expansions are:<br>
<br>
 asin (e)     -> e<br>
 sqrt (1-e^2) -> 1 - e^2/2<br>
<br>
When you insert these limits into the AA3 equation, then:<br>
<br>
 AA3 -> 2/e^2 - 2*(1-e^2/2)*e / e^3<br>
     -> 2/e^2 - 2/e^2 + 2*e^2/2<br>
     -> 1<br>
<br>
so there is no error in the formula.<br>
<br>
Cheers,<br>
Norbert<br>
<br>
<br>
On Wed, 14 Aug 2013, Chenchong Zhu wrote:<br>
<br>
 Hi FLASH users,<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
I'm currently playing with the MacLaurin spheroid example that comes with<br>
FLASH to see how well the spheroid can be represented hydrodynamically.<br>
I'm not entirely certain, but I think I might have found two bugs in the<br>
code though.<br>
<br>
On line 78 of Simulation_init.F90 AA3 as written approaches 2/3 rather<br>
than<br>
1 as sim_eccentricity goes to 0.  It obeys the proper limit if it's<br>
changed<br>
to "AA3 = (3.0*sqrt...".<br>
<br>
In Simulation_initBlock.F90, under case(sim_geom3DCartesian), vxfac and<br>
vyfac have z-dependencies on their velocities.  They should be multiplied<br>
by rxyinv = 1./sqrt(xdist**2 + ydist**2) rather than rinv if we want<br>
rigid<br>
rotation.<br>
<br>
Best,<br>
<br>
Charles<br>
<br>
<br>
</blockquote>
<br>
</blockquote></blockquote>
<br>
</blockquote>
</div></div></blockquote></div><br></div>