[FLASH-USERS] Minor bugs in MacLaurin spheroid test

[Norbert Flocke]

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
>>>>
>>>>
>>>
>