[CIG-MC] Re: geoid scaling in CitcomS-2.2.1
Eh Tan
tan2 at geodynamics.org
Mon Oct 29 14:59:08 PDT 2007
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
More information about the CIG-MC
mailing list