[cig-commits] r5956 - in mc/3D/CitcomS/trunk: CitcomS/Components
lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Feb 2 15:30:07 PST 2007
Author: tan2
Date: 2007-02-02 15:30:07 -0800 (Fri, 02 Feb 2007)
New Revision: 5956
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/IC.py
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Initial_temperature.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
A new option for initial temperature, which is used in Shijie's benchmark
When tic_method==3, the initial temperature is a linear profile plus
mag * modified_plgndr_a(ll, mm, theta) * (cos(mm*phi) + sin(mm*phi))
* sin(M_PI * (r - E->sphere.ri) / (E->sphere.ro - E->sphere.ri))
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/IC.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/IC.py 2007-02-02 22:54:56 UTC (rev 5955)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/IC.py 2007-02-02 23:30:07 UTC (rev 5956)
@@ -106,9 +106,9 @@
zero_elapsed_time = pyre.inventory.bool("zero_elapsed_time", default=True)
tic_method = pyre.inventory.int("tic_method", default=0,
- validator=pyre.inventory.choice([0, 1, 2]))
+ validator=pyre.inventory.choice([0, 1, 2, 3]))
- # for tic_method=0
+ # for tic_method=0 or 3
num_perturbations = pyre.inventory.int("num_perturbations", default=1,
validator=pyre.inventory.less(255))
perturbl = pyre.inventory.list("perturbl", default=[1])
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-02-02 22:54:56 UTC (rev 5955)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-02-02 23:30:07 UTC (rev 5956)
@@ -308,6 +308,46 @@
}
}
}
+ else if (E->convection.tic_method == 3) {
+ /* set up a linear temperature profile first */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz;
+ r1=E->sx[m][3][node];
+ E->T[m][node] = E->control.TBCbotval - (E->control.TBCtopval + E->control.TBCbotval)*(r1 - E->sphere.ri)/(E->sphere.ro - E->sphere.ri);
+ }
+
+ /* This part put a temperature anomaly at depth where the global
+ node number is equal to load_depth. The horizontal pattern of
+ the anomaly is given by spherical harmonic ll & mm. */
+
+ for (p=0; p<E->convection.number_of_perturbations; p++) {
+ mm = E->convection.perturb_mm[p];
+ ll = E->convection.perturb_ll[p];
+ con = E->convection.perturb_mag[p];
+ kk = E->convection.load_depth[p];
+
+ if ( (kk < 1) || (kk >= gnoz) ) continue;
+
+ if (E->parallel.me == 0)
+ fprintf(stderr,"Initial temperature perturbation: layer=%d mag=%g l=%d m=%d\n", kk, con, ll, mm);
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz;
+ t1=E->sx[m][1][node];
+ f1=E->sx[m][2][node];
+ r1=E->sx[m][3][node];
+ E->T[m][node] += con*modified_plgndr_a(ll,mm,t1)
+ *(cos(mm*f1)+sin(mm*f1))
+ *sin(M_PI*(r1-E->sphere.ri)/(E->sphere.ro-E->sphere.ri));
+ }
+ }
+ }
else if (E->convection.tic_method == 1) {
}
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-02-02 22:54:56 UTC (rev 5955)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2007-02-02 23:30:07 UTC (rev 5956)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
@@ -64,10 +64,14 @@
- blob_center: A comma-separated list of three float numbers.
- blob_radius: A dmensionless length, typically a fraction
of the Earth's radius.
- - blob_dT : Dimensionless temperature. */
+ - blob_dT : Dimensionless temperature.
+ When tic_method is 3, the temperature is a linear profile + perturbation
+ for whole mantle.
+ */
- if (E->convection.tic_method == 0) {
+
+ if (E->convection.tic_method == 0 || E->convection.tic_method == 3 ) {
/* This part put a temperature anomaly at depth where the global
node number is equal to load_depth. The horizontal pattern of
the anomaly is given by spherical harmonic ll & mm. */
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-02-02 22:54:56 UTC (rev 5955)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-02-02 23:30:07 UTC (rev 5956)
@@ -445,7 +445,45 @@
E->T[m][node] += amp * exp(-1.0*distance/radius);
}
}
+ else if (E->convection.tic_method == 3) {
+ /* set up a linear temperature profile first */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz;
+ r1=E->sx[m][3][node];
+ E->T[m][node] = E->control.TBCbotval - (E->control.TBCtopval + E->control.TBCbotval)*(r1 - E->sphere.ri)/(E->sphere.ro - E->sphere.ri);
+ }
+ /* This part put a temperature anomaly for whole mantle. The horizontal
+ pattern of the anomaly is given by spherical harmonic ll & mm. */
+
+ for (p=0; p<E->convection.number_of_perturbations; p++) {
+ mm = E->convection.perturb_mm[p];
+ ll = E->convection.perturb_ll[p];
+ con = E->convection.perturb_mag[p];
+ kk = E->convection.load_depth[p];
+
+ if ( (kk < 1) || (kk >= gnoz) ) continue;
+
+ if (E->parallel.me == 0)
+ fprintf(stderr,"Initial temperature perturbation: layer=%d mag=%g l=%d m=%d\n", kk, con, ll, mm);
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz;
+ t1=E->sx[m][1][node];
+ f1=E->sx[m][2][node];
+ r1=E->sx[m][3][node];
+ E->T[m][node] += con*(cos(mm*f1)+sin(mm*f1))
+ *sin(M_PI*(r1-E->sphere.ri)/(E->sphere.ro-E->sphere.ri));
+ }
+ }
+ }
+
temperatures_conform_bcs(E);
return;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-02-02 22:54:56 UTC (rev 5955)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-02-02 23:30:07 UTC (rev 5956)
@@ -237,7 +237,7 @@
getIntProperty(properties, "tic_method", E->convection.tic_method, fp);
- if (E->convection.tic_method == 0) {
+ if (E->convection.tic_method == 0 || E->convection.tic_method == 3) {
int num_perturb;
getIntProperty(properties, "num_perturbations", num_perturb, fp);
More information about the cig-commits
mailing list