[FLASH-USERS] Minor bugs in MacLaurin spheroid test

[Chenchong Zhu]

Hi Norbert,

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:

  p0     = (2.0/3.0)*sim_pi*sim_Newton*sim_density**2*amean**2
  sim_Pconst = p0 * AA3 * (1.0-sim_eccentricity**2)**(2.0/3.0)


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.

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.

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?

Thanks again for your help,

Charles


On Thu, Aug 15, 2013 at 11:07 AM, Norbert Flocke
<flocke at flash.uchicago.edu>wrote:

> Hi Charles,
>
> You are right, I did not consider enough terms. On the other hand, the
> AA3 going to 2/3 with e -> 0 is the correct behavior. Please look at the
> 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
> two constants A1 and A3. When e -> 0 both go to 2/3 and the result is the
> potential of a uniform density sphere.
>
> 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/**Multipole).
> Please have a look at the routine:
>
>    sim_initBlockAnalytical.F90
>
> 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.
>
>
> Cheers,
> Norbert
>
>
>
>
> On Wed, 14 Aug 2013, Chenchong Zhu wrote:
>
>  Hi Norbert,
>>
>> Thanks for the response!
>>
>> For sqrt(1-e^2)*asin (e), I don't think you kept enough of the expansion:
>>
>> asin(e) -> e + e^3/6
>> sqrt(1 - e^2) -> 1 - e^2/2 - e^4/8 -> 1 - e^2/2
>>
>> AA3 = 2/e^2 - 2*sqrt(1-e^2)*asin (e) / e^3
>>       = 2/e^2 - (2/e^3)*(e - e^3/3 - e^5/12)
>>       = 2/e^2 - 2/e^2 + 2/3
>>       = 2/3
>>
>> Throwing the expression into Python and setting e = tiny number gives the
>> same result.
>>
>> Best,
>>
>> Charles
>>
>>
>> On Wed, Aug 14, 2013 at 1:29 PM, Norbert Flocke
>> <flocke at flash.uchicago.edu>**wrote:
>>
>>  Sorry, there is a typo in the last AA3 equations. It should read:
>>>
>>>
>>>   AA3 -> 2/e^2 - 2*(1-e^2/2)*e / e^3
>>>       -> 2/e^2 - 2/e^2 + (2*e^2/2*e^2)
>>>       -> 1
>>>
>>> Norbert
>>>
>>>
>>> On Wed, 14 Aug 2013, Norbert Flocke wrote:
>>>
>>>  Hi Charles,
>>>
>>>>
>>>> The answer to your first question I can provide. As written in the code
>>>> on line 78, the value of AA3 is
>>>>
>>>>  AA3 = (2*sqrt(1-e^2)/e^2) * (1/sqrt(1-e^2) - asin (e)/e)
>>>>
>>>> or, written out
>>>>
>>>>  AA3 = 2/e^2 - 2*sqrt(1-e^2)*asin (e) / e^3
>>>>
>>>> where e is the eccentricity. If e -> zero, then the MacLaurin (what
>>>> coincidence!) expansions are:
>>>>
>>>>  asin (e)     -> e
>>>>  sqrt (1-e^2) -> 1 - e^2/2
>>>>
>>>> When you insert these limits into the AA3 equation, then:
>>>>
>>>>  AA3 -> 2/e^2 - 2*(1-e^2/2)*e / e^3
>>>>      -> 2/e^2 - 2/e^2 + 2*e^2/2
>>>>      -> 1
>>>>
>>>> so there is no error in the formula.
>>>>
>>>> Cheers,
>>>> Norbert
>>>>
>>>>
>>>> On Wed, 14 Aug 2013, Chenchong Zhu wrote:
>>>>
>>>>  Hi FLASH users,
>>>>
>>>>>
>>>>> I'm currently playing with the MacLaurin spheroid example that comes
>>>>> with
>>>>> FLASH to see how well the spheroid can be represented hydrodynamically.
>>>>> I'm not entirely certain, but I think I might have found two bugs in
>>>>> the
>>>>> code though.
>>>>>
>>>>> On line 78 of Simulation_init.F90 AA3 as written approaches 2/3 rather
>>>>> than
>>>>> 1 as sim_eccentricity goes to 0.  It obeys the proper limit if it's
>>>>> changed
>>>>> to "AA3 = (3.0*sqrt...".
>>>>>
>>>>> In Simulation_initBlock.F90, under case(sim_geom3DCartesian), vxfac and
>>>>> vyfac have z-dependencies on their velocities.  They should be
>>>>> multiplied
>>>>> by rxyinv = 1./sqrt(xdist**2 + ydist**2) rather than rinv if we want
>>>>> rigid
>>>>> rotation.
>>>>>
>>>>> Best,
>>>>>
>>>>> Charles
>>>>>
>>>>>
>>>>>
>>>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20130815/fd5034c8/attachment.htm>