[CIG-MC] Re: geoid scaling in CitcomS-2.2.1

[Eh Tan]

Scott King wrote:
> I've been checking the scaling in CitcomS.
>
> I think it is o.k. but the the way it is done is a bit confusing.
>
> Look in Topo_Gravity.c
>
> in geoid_from_buoyancy:
>
>      /* scale for buoyancy */
>     scaling2 =
> -E->data.therm_exp*E->data.ref_temperature*E->data.density/E->control.Atemp;
>
>     /* scale for geoid */
>     scaling = 1.0e6*4.0*M_PI*E->data.radius_km*E->data.radius_km*
>         E->data.grav_const/E->data.grav_acc;
>
> in geoid_from_topography:
>
>     stress_scaling = E->data.ref_viscosity*E->data.therm_diff/
>         (E->data.radius_km*E->data.radius_km*1e6);
>
>     /* density contrast across surface */
>     den_contrast1 = (E->data.density-E->data.density_above);
>     /* density contrast across CMB */
>     den_contrast2 = (E->data.density_below-E->data.density);
>
>     /* scale for surface and CMB topo */
>     topo_scaling1 = stress_scaling / (den_contrast1*E->data.grav_acc);
>     topo_scaling2 = stress_scaling / (den_contrast2*E->data.grav_acc);
>
>     /* scale for geoid */
>     scaling = 1.0e3*4.0*M_PI*E->data.radius_km*
>         E->data.grav_const/E->data.grav_acc;
>
>
> I would suggest that it would be more clear if in both routines, the
> "scaling for geoid term" term  looked like:
>
>     /* scale for geoid */
>     scaling = 1.0e3*4.0*M_PI*E->data.radius_km*
>         E->data.grav_const/E->data.grav_acc;
>
> This would mean in geoid_from_buoyancy
>     scaling2 =
> -1.0e3*E->data.therm_exp*E->data.ref_temperature*E->data.density*E->data.radius_km/E->control.Atemp;
>
>
> It took me a while to figure out what was going on.  I did not notice
> that scaling was different in the two terms and it looked
> like there was an error.

Hi Scott,

(Cc'ing to cig-mc as well.)

You are right. The "scale for geoid" in those two functions are referred
to different constants. Using the same comment there is quite
misleading. I will clarify the code and provide the actual equations in
the comment soon.


>
> I still don't understand the strange value I'm getting in the
> geoid_from_topography for one term (all other terms are o.k. except
> for l=m=2
> 2 2 1.8237e+28 1.3652e-03 -4.4414e+01 1.2344e-03 1.8237e+28 1.5026e-04
>
> For other problems (same grid) the problem seems to be in
> l=3,m=1...    It is always in the stress term of the geoid.  It is like
> there is something nasty with the stresses.    Do you ever play with
> the parameters for the convergence when you are looking
> at stresses?   Something is odd. 

This sounds strange. I haven't encounter this problem before.

My guess is that your velocity solution somehow goes haywire. Sometimes,
the Stokes' solver has a hard time to converge, and the symptoms are
surprisingly large amplitude in velocity (>10^3 larger than it should
be) and a large iteration number for Uzawa algorithm. This usually
occurs when the viscosity changes rapidly in space or when complex plate
motion is imposed as boundary conditions. Is this the case?


-- 
Eh Tan
Staff Scientist
Computational Infrastructure for Geodynamics
2750 E. Washington Blvd. Suite 210
Pasadena, CA 91107
(626) 395-1693
http://www.geodynamics.org