[CIG-MC] Ridge rots and blow ups in CitcomS
Thorsten Becker
twb at usc.edu
Sat Oct 13 08:37:43 PDT 2007
Dear Shijie,
Thanks for the explanation, great to hear that there's a benchmark
going on. I am working on the SVN version of CitcomS when adding my
modifications, that's why I got the update already.
As you know, following your 2001 EPSL paper, I'm a big fan of the net
rotations myself, and noticed that those can be quite strong
particularly for plastic yielding spherical models, quite interesting.
That's why I put in the lib/Determine_net_rotation.c a while back.
However, as to be expected, your implementation of the net rotation
removal should be much quicker than mine. Also, I forgot about the time
step issue and agree that it's much better to remove the net rotation
after every velocity solution, rather than removing only at output.
The problem that caused the blowup when your routine was called within
CIG CitcomS was very simple; I would fix it in the SVN version but
cannot access my cluster right now:
VV should be declared float rather than double in remove_rigid_rot
float VV[4][9]; (l. 819 of Global_operations.c)
because (at least in my version)
void velo_from_element(E,VV,m,el,sphere_key)
struct All_variables *E;
float VV[4][9];
int el,m,sphere_key;
in Nodal_mesh.c
It would be very easy to allow the compiler to catch those declarative
mismatches from now on for eternity, if we would switch the function
declarations from old-school type to the
void velo_from_element(struct All_variables *,
float [4][9],int ,int ,int);
type. Each file could include a "functions.h" file which has all the
function declarations. This could be generated automatically by
cproto $(INCLUDES) -f2 -q *.c | grep -v "void main(" | grep -v "int
main(" > functions.h
In fact, Luiz once made such changes, but then they must have been taken
out of the major source tree of CitcomS for some reason.
Cheers
T
On Sat, 2007-10-13 at 03:18 -0600, Shijie Zhong wrote:
> Thorsten,
>
> remove_rigid_rot() was added only ~ a week ago (you really keep things updated). In the new release (possibly next week?), the purpose of this routine will be explained. I am working on a benchmark paper (would have been done already if not for this rigid body rotation) that would detail this issue. I can explain to you (and those who may be interested) briefly here.
>
> For the stokes' flow in a spherical shell, it is possible to add an arbitrary ridge body rotation to the flow and the governing equations will still be satisfied (i.e., the rigid body rotation is unconstrained). In most calculations (as far as I know), this mode of rigid body rotation is never important. However, for some reasons, it may show up in some other cases (e.g., thermochemical convection, as Allen McNamara discovered, or mobile-lid convection with continental keels). When the rigid body rotation is present in a calculation, it can seriously and unnecessarily reduce time step (delta t) in convection calculation, as delta t is determined by maximum velocity, and it also leads to ambiguity in velocities (i.e., reference frame issue, important for plate tectonics related topics and benchmark purposes).
>
> I think that a proper reference frame is no-net-rotation for the whole MANTLE (different from no-net-rotation for lithosphere as used sometimes in plate tectonics related studies). In the past, as what you have done, I have used a post-processing routine to remove the rigid body rotation for the whole mantle (e.g., Zhong [2001]). Allen then found that the rigid body rotation sometimes gave him too small delta t in time-dependent thermochemical convection calculations, so he coded a routine into CitcomS (thermochemical one in 2004) that can be called every time step.
>
> I recently "re-learnt" the lesson when working on a benchmark paper (more complete description on this issue is in that paper). I "wasted" more than a month in trying to figure out why two resolution test cases for thermochemical convection that differ only in resolution have very different Vrms, while heat flux is nearly the same. It turns out that the high resolution case somehow has some rigid body rotation.
> At the end, I think that we should include this routine and that the rigid body rotation should be removed every time step, although it does not affect the dynamics. This also gives us a reference frame for all the velocities (i.e., no-net-rotation for the mantle). There should be two routines in CitcomS right now, one is Allen's original one, and the other (being called) is a simplified version from Allen's.
>
> Now to the non-Newtonian issue, I cannot see why this routine should not mess up the non-Newtonian iterations. There might be some incompartibility issue between this routine and CIG's version (?). Perhaps, you can turn it off for now.
>
>
> Shijie
>
>
> Shijie Zhong
> Department of Physics
> University of Colorado at Boulder
> Boulder, CO 80309
> Tel: 303-735-5095; Fax: 303-492-7935
> Web: http://anquetil.colorado.edu/szhong
>
>
> ---- Original message ----
> >Date: Fri, 12 Oct 2007 17:59:15 -0700
> >From: Thorsten Becker <twb at usc.edu>
> >Subject: [CIG-MC] Ridge rots and blow ups in CitcomS
> >To: Eh Tan <tan2 at geodynamics.org>
> >Cc: cig-mc at geodynamics.org
> >
> >Hi Eh,
> >
> >I am a bit confused about the net rotation removal routine
> >
> >remove_rigid_rot in Global_operations.c
> >
> >in CitcomS. Has this always been there, or is this a relatively new
> >addition?
> >
> >I had recently also built in a (slightly different) net rotation
> >removal routine as
> >
> >lib/Determine_net_rotation.c:53:double
> >determine_model_net_rotation(struct All_variables *,double *);
> >
> >which I activate only in the gzipped output routine in a post-processing
> >step if
> >
> >gzdir_rnr=on
> >
> >
> >and not at every time step. I noticed the new routine because it lead
> >to a blow up of velocities if stress-dependent viscosity iterations
> >were performed after a tic_method=-1 restart, for reasons unclear to me
> >right now.
> >
> >Could we maybe set
> >
> >remove_rigid_rotation
> >
> >to an "off" default for now?
> >
> >Thanks for any comments
> >
> >T
> >
> >
> >
> >
> >--
> > Thorsten W Becker Department of Earth Sciences
> > University of Southern California p: 213.740.8365
> > Los Angeles CA 90089-0740 http://geodynamics.usc.edu/~becker
> >
> >_______________________________________________
> >CIG-MC mailing list
> >CIG-MC at geodynamics.org
> >http://geodynamics.org/cgi-bin/mailman/listinfo/cig-mc
--
Thorsten W Becker Department of Earth Sciences
University of Southern California p: 213.740.8365
Los Angeles CA 90089-0740 http://geodynamics.usc.edu/~becker
More information about the CIG-MC
mailing list