[FLASH-USERS] Minor bugs in MacLaurin spheroid test

[Chenchong Zhu]

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/20130814/d0119ba5/attachment.htm>