[cig-commits] r13273 - in mc/3D/CitcomS/trunk: CitcomS/Components examples/Cookbook6 lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Nov 7 15:32:20 PST 2008
Author: tan2
Date: 2008-11-07 15:32:20 -0800 (Fri, 07 Nov 2008)
New Revision: 13273
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/IC.py
mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg
mc/3D/CitcomS/trunk/lib/Full_solver.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Initial_temperature.c
mc/3D/CitcomS/trunk/lib/Regional_solver.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/solver.h
Log:
* Reserved 'tic_method=100' for user-defined initial temperature.
* Read all parameters related to initial temperature regardless which tic_method is used. This will make the life easier when adding new tic_method.
* Fixed a bug in 'tic_method=0' for regional model. The bug causes that the sinosoidal temperature perturbation is applied at 0th processor only.
* Fixed a bug in 'tic_method=1' for regional model. 'tic_method=1' should generate a top TBL according to the input parameter 'half_space_age'. However, the code effectively multiply 25x to the age. Cookbook6 is using this tic_method, so its half_space_age needs to become 2500 to get the original temperature.
* In 'tic_method=2', compute distance in Cartesian coordinate, instead of spherical coordinate.
* Refactoring codes for temperature initial conditions.
* Merged regional_construct_tic_from_input() and full_construct_tic_from_input() to construct_tic_from_input().
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/IC.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/IC.py 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/IC.py 2008-11-07 23:32:20 UTC (rev 13273)
@@ -115,8 +115,7 @@
solution_cycles_init = pyre.inventory.int("solution_cycles_init", default=0)
zero_elapsed_time = pyre.inventory.bool("zero_elapsed_time", default=True)
- tic_method = pyre.inventory.int("tic_method", default=0,
- validator=pyre.inventory.choice([-1, 0, 1, 2, 3]))
+ tic_method = pyre.inventory.int("tic_method", default=0)
# for tic_method=0 or 3
num_perturbations = pyre.inventory.int("num_perturbations", default=1,
Modified: mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg 2008-11-07 23:32:20 UTC (rev 13273)
@@ -48,7 +48,7 @@
blob_center = 1.570800e+00,1.570800e+00,9.246600e-01
blob_radius = 6.278334e-02
blob_dT = 0.18
-half_space_age = 100.0
+half_space_age = 2500.0
[CitcomS.solver.bc]
Modified: mc/3D/CitcomS/trunk/lib/Full_solver.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_solver.c 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/lib/Full_solver.c 2008-11-07 23:32:20 UTC (rev 13273)
@@ -56,7 +56,6 @@
/* Version_dependent.c */
void full_node_locations(struct All_variables *);
-void full_construct_tic_from_input(struct All_variables *);
void full_construct_boundary(struct All_variables *);
@@ -88,7 +87,6 @@
/* Version_dependent.c */
E->solver.node_locations = full_node_locations;
- E->solver.construct_tic_from_input = full_construct_tic_from_input;
E->solver.construct_boundary = full_construct_boundary;
return;
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2008-11-07 23:32:20 UTC (rev 13273)
@@ -232,137 +232,6 @@
-void full_construct_tic_from_input(struct All_variables *E)
-{
- int i, j, k, kk, m, p, node;
- int nox, noy, noz, gnoz;
- double r1, f1, t1;
- int mm, ll;
- double con;
- double modified_plgndr_a(int, int, double);
- void debug_sphere_expansion(struct All_variables *);
- void temperatures_conform_bcs();
-
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- gnoz=E->mesh.noz;
-
-
- switch (E->convection.tic_method){
-
- case 0:
-
- /* 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;
-
- k = kk - E->lmesh.nzs + 1;
- if ( (k < 1) || (k >= noz) ) continue; /* if layer k is not inside this proc. */
- if (E->parallel.me_loc[1] == 0 && E->parallel.me_loc[2] == 0
- && E->sphere.capid[1] == 1 )
- 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++) {
- node=k+(j-1)*noz+(i-1)*nox*noz;
- t1=E->sx[m][1][node];
- f1=E->sx[m][2][node];
-
- E->T[m][node] += con*modified_plgndr_a(ll,mm,t1)*cos(mm*f1);
- }
- }
- break;
-
- case 3:
-
- /* set up a conductive 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.TBCtopval*E->sphere.ro
- - E->control.TBCbotval*E->sphere.ri)
- / (E->sphere.ro - E->sphere.ri)
- + (E->control.TBCbotval - E->control.TBCtopval)
- * E->sphere.ro * E->sphere.ri / r1
- / (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*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));
- }
- }
- break;
- case -1:
- case 1:
- case 2:
- break;
- case 4: /* read initial temp from grd grd files */
-#ifdef USE_GGRD
- ggrd_full_temp_init(E);
-#else
- fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
- parallel_process_termination();
-#endif
- break;
- default: /* unknown option */
- fprintf(stderr,"Invalid value of 'tic_method'\n");
- parallel_process_termination();
- break;
- }
-
- temperatures_conform_bcs(E);
-
- /* debugging the code of expanding spherical harmonics */
- /* debug_sphere_expansion(E);*/
- return;
-}
-
-
/* setup boundary node and element arrays for bookkeeping */
void full_construct_boundary( struct All_variables *E)
Modified: mc/3D/CitcomS/trunk/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/lib/Initial_temperature.c 2008-11-07 23:32:20 UTC (rev 13273)
@@ -28,17 +28,20 @@
#include <math.h>
#include <assert.h>
+#include <string.h>
#include "global_defs.h"
#include "lith_age.h"
#include "parsing.h"
void parallel_process_termination();
-void temperatures_conform_bcs();
+void temperatures_conform_bcs(struct All_variables *);
+double modified_plgndr_a(int, int, double);
#include "initial_temperature.h"
-void debug_tic(struct All_variables *);
-void read_tic_from_file(struct All_variables *);
+static void debug_tic(struct All_variables *);
+static void read_tic_from_file(struct All_variables *);
+static void construct_tic_from_input(struct All_variables *);
#ifdef USE_GZDIR
void restart_tic_from_gzdir_file(struct All_variables *);
@@ -53,7 +56,7 @@
int m = E->parallel.me;
int noz = E->lmesh.noz;
- int n,tmp;
+ int n;
input_int("tic_method", &(E->convection.tic_method), "0,0,2", m);
@@ -90,11 +93,6 @@
*/
- switch(E->convection.tic_method){
- case -1: /* read from file, no other options needed */
- break;
- case 0:
- case 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. */
@@ -129,14 +127,8 @@
E->convection.load_depth[0] = (noz+1)/2;
}
- break;
- case 1: /* case 1 */
-
input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
- break;
- case 2: /* case 2 */
- input_float("half_space_age", &(E->convection.half_space_age), "40.0,1e-3,nomax", m);
if( ! input_float_vector("blob_center", 3, E->convection.blob_center, m)) {
assert( E->sphere.caps == 12 || E->sphere.caps == 1 );
if(E->sphere.caps == 12) { /* Full version: just quit here */
@@ -152,8 +144,7 @@
}
input_float("blob_radius", &(E->convection.blob_radius), "0.063,0.0,1.0", m);
input_float("blob_dT", &(E->convection.blob_dT), "0.18,nomin,nomax", m);
- break;
- case 4:
+
/*
case 4: initial temp from grd files
*/
@@ -181,13 +172,6 @@
parallel_process_termination();
#endif
- break;
- default: /* unknown option */
- fprintf(stderr,"Invalid value of 'tic_method'\n");
- parallel_process_termination();
- break;
- }
-
return;
}
@@ -209,9 +193,9 @@
read_tic_from_file(E);
}
else if (E->control.lith_age)
- lith_age_construct_tic(E);
+ lith_age_construct_tic(E);
else
- (E->solver.construct_tic_from_input)(E);
+ construct_tic_from_input(E);
/* Note: it is the callee's responsibility to conform tbc. */
/* like a call to temperatures_conform_bcs(E); */
@@ -223,7 +207,7 @@
}
-void debug_tic(struct All_variables *E)
+static void debug_tic(struct All_variables *E)
{
int m, j;
@@ -240,10 +224,8 @@
-void read_tic_from_file(struct All_variables *E)
+static void read_tic_from_file(struct All_variables *E)
{
- void temperatures_conform_bcs();
-
int ii, ll, mm;
float tt;
int i, m;
@@ -286,3 +268,404 @@
}
+static void linear_temperature_profile(struct All_variables *E)
+{
+ int m, i, j, k, node;
+ int nox, noy, noz;
+ double r1;
+
+ nox = E->lmesh.nox;
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+
+ 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);
+ }
+
+ return;
+}
+
+
+static void conductive_temperature_profile(struct All_variables *E)
+{
+ int m, i, j, k, node;
+ int nox, noy, noz;
+ double r1;
+
+ nox = E->lmesh.nox;
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+
+ 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.TBCtopval*E->sphere.ro
+ - E->control.TBCbotval*E->sphere.ri)
+ / (E->sphere.ro - E->sphere.ri)
+ + (E->control.TBCbotval - E->control.TBCtopval)
+ * E->sphere.ro * E->sphere.ri / r1
+ / (E->sphere.ro - E->sphere.ri);
+ }
+
+ return;
+}
+
+
+static void constant_temperature_profile(struct All_variables *E, double mantle_temp)
+{
+ int m, i;
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=1; i<=E->lmesh.nno; i++)
+ E->T[m][i] = mantle_temp;
+
+ return;
+}
+
+
+static void add_top_tbl(struct All_variables *E, double age_in_myrs, double mantle_temp)
+{
+ int m, i, j, k, node;
+ int nox, noy, noz;
+ double r1, dT, tmp;
+
+ nox = E->lmesh.nox;
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+
+ dT = (mantle_temp - E->control.TBCtopval);
+ tmp = 0.5 / sqrt(age_in_myrs / E->data.scalet);
+
+ fprintf(stderr, "%e %e\n", dT, tmp);
+ 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] -= dT * erfc(tmp * (E->sphere.ro - r1));
+ }
+
+ return;
+}
+
+
+static void add_bottom_tbl(struct All_variables *E, double age_in_myrs, double mantle_temp)
+{
+ int m, i, j, k, node;
+ int nox, noy, noz;
+ double r1, dT, tmp;
+
+ nox = E->lmesh.nox;
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+
+ dT = (E->control.TBCbotval - mantle_temp);
+ tmp = 0.5 / sqrt(age_in_myrs / E->data.scalet);
+
+ 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] += dT * erfc(tmp * (r1 - E->sphere.ri));
+ }
+
+ return;
+}
+
+
+static void add_perturbations_at_layers(struct All_variables *E)
+{
+ /* This function 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 wavenumber (in regional model) or
+ by spherical harmonic (in global model). */
+
+ int m, i, j, k, node;
+ int p, ll, mm, kk;
+ int nox, noy, noz, gnoz;
+ double t1, f1, tlen, flen, con;
+
+ nox = E->lmesh.nox;
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+ gnoz = E->mesh.noz;
+
+ for (p=0; p<E->convection.number_of_perturbations; p++) {
+ ll = E->convection.perturb_ll[p];
+ mm = E->convection.perturb_mm[p];
+ kk = E->convection.load_depth[p];
+ con = E->convection.perturb_mag[p];
+
+ if ( (kk < 1) || (kk >= gnoz) ) continue; /* layer kk is outside domain */
+
+ k = kk - E->lmesh.nzs + 1; /* convert global nz to local nz */
+ if ( (k < 1) || (k >= noz) ) continue; /* layer k is not inside this proc. */
+ if (E->parallel.me_loc[1] == 0 && E->parallel.me_loc[2] == 0
+ && E->sphere.capid[1] == 1 )
+ fprintf(stderr,"Initial temperature perturbation: layer=%d mag=%g l=%d m=%d\n", kk, con, ll, mm);
+
+ if(E->sphere.caps == 1) {
+ /* regional mode, add sinosoidal perturbation */
+
+ tlen = M_PI / (E->control.theta_max - E->control.theta_min);
+ flen = M_PI / (E->control.fi_max - E->control.fi_min);
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=1; i<=noy; i++)
+ for(j=1; j<=nox;j ++) {
+ node = k + (j-1)*noz + (i-1)*nox*noz;
+ t1 = (E->sx[m][1][node] - E->control.theta_min) * tlen;
+ f1 = (E->sx[m][2][node] - E->control.fi_min) * flen;
+
+ E->T[m][node] += con * cos(ll*t1) * cos(mm*f1);
+ }
+ }
+ else {
+ /* global mode, add spherical harmonics perturbation */
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=1; i<=noy; i++)
+ for(j=1; j<=nox;j ++) {
+ node = k + (j-1)*noz + (i-1)*nox*noz;
+ t1 = E->sx[m][1][node];
+ f1 = E->sx[m][2][node];
+
+ E->T[m][node] += con * modified_plgndr_a(ll,mm,t1) * cos(mm*f1);
+ }
+ } /* end if */
+ } /* end for p */
+
+ return;
+}
+
+
+static void add_perturbations_at_all_layers(struct All_variables *E)
+{
+ /* This function put a temperature anomaly for whole mantle with
+ a sinosoidal amplitude in radial dependence. The horizontal pattern
+ of the anomaly is given by wavenumber (in regional model) or
+ by spherical harmonic (in global model). */
+
+ int m, i, j, k, node;
+ int p, ll, mm;
+ int nox, noy, noz, gnoz;
+ double r1, t1, f1, tlen, flen, rlen, con;
+
+ nox = E->lmesh.nox;
+ noy = E->lmesh.noy;
+ noz = E->lmesh.noz;
+ gnoz = E->mesh.noz;
+
+ rlen = M_PI / (E->sphere.ro - E->sphere.ri);
+
+ for (p=0; p<E->convection.number_of_perturbations; p++) {
+ ll = E->convection.perturb_ll[p];
+ mm = E->convection.perturb_mm[p];
+ con = E->convection.perturb_mag[p];
+
+ if (E->parallel.me_loc[1] == 0 && E->parallel.me_loc[2] == 0
+ && E->sphere.capid[1] == 1 )
+ fprintf(stderr,"Initial temperature perturbation: mag=%g l=%d m=%d\n", con, ll, mm);
+
+ if(E->sphere.caps == 1) {
+ /* regional mode, add sinosoidal perturbation */
+
+ tlen = M_PI / (E->control.theta_max - E->control.theta_min);
+ flen = M_PI / (E->control.fi_max - E->control.fi_min);
+
+ 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] - E->control.theta_min) * tlen;
+ f1 = (E->sx[m][2][node] - E->control.fi_min) * flen;
+ r1 = E->sx[m][3][node];
+
+ E->T[m][node] += con * cos(ll*t1) * cos(mm*f1)
+ * sin((r1-E->sphere.ri) * rlen);
+ }
+ }
+ else {
+ /* global mode, add spherical harmonics perturbation */
+
+ 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((r1-E->sphere.ri) * rlen);
+ }
+ } /* end if */
+ } /* end for p */
+
+ return;
+}
+
+
+static void add_spherical_anomaly(struct All_variables *E)
+{
+ int i, j ,k , m, node;
+ int nox, noy, noz;
+
+ double theta_center, fi_center, r_center;
+ double radius, amp;
+
+ double x_center, y_center, z_center;
+ double x, y, z, distance;
+
+ noy = E->lmesh.noy;
+ nox = E->lmesh.nox;
+ noz = E->lmesh.noz;
+
+ theta_center = E->convection.blob_center[0];
+ fi_center = E->convection.blob_center[1];
+ r_center = E->convection.blob_center[2];
+ radius = E->convection.blob_radius;
+ amp = E->convection.blob_dT;
+
+ fprintf(stderr,"center=(%e %e %e) radius=%e dT=%e\n",
+ theta_center, fi_center, r_center, radius, amp);
+
+ x_center = r_center * sin(fi_center) * cos(theta_center);
+ y_center = r_center * sin(fi_center) * sin(theta_center);
+ z_center = r_center * cos(fi_center);
+
+ /* compute temperature field according to nodal coordinate */
+ 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;
+
+ x = E->x[m][1][node];
+ y = E->x[m][2][node];
+ z = E->x[m][3][node];
+
+ distance = sqrt((x-x_center)*(x-x_center) +
+ (y-y_center)*(y-y_center) +
+ (z-z_center)*(z-z_center));
+
+ if (distance < radius)
+ E->T[m][node] += amp * exp(-1.0*distance/radius);
+ }
+ return;
+}
+
+
+static void construct_tic_from_input(struct All_variables *E)
+{
+ double mantle_temperature;
+
+ switch (E->convection.tic_method){
+ case 0:
+ /* a linear temperature profile + perturbations at some layers */
+ linear_temperature_profile(E);
+ add_perturbations_at_layers(E);
+ break;
+
+ case 1:
+ /* T=1 for whole mantle + cold lithosphere TBL */
+ mantle_temperature = 1;
+ constant_temperature_profile(E, mantle_temperature);
+ add_top_tbl(E, E->convection.half_space_age, mantle_temperature);
+ break;
+
+ case 2:
+ /* T='mantle_temp' for whole mantle + cold lithosphere TBL
+ + a spherical anomaly at lower center */
+ mantle_temperature = E->control.lith_age_mantle_temp;
+ constant_temperature_profile(E, mantle_temperature);
+ add_top_tbl(E, E->convection.half_space_age, mantle_temperature);
+ add_spherical_anomaly(E);
+ break;
+
+ case 3:
+ /* a conductive temperature profile + perturbations at all layers */
+ conductive_temperature_profile(E);
+ add_perturbations_at_all_layers(E);
+ break;
+
+ case 4:
+ /* read initial temperature from grd files */
+#ifdef USE_GGRD
+ if (E->sphere.caps == 1)
+ ggrd_reg_temp_init(E);
+ else
+ ggrd_full_temp_init(E);
+#else
+ fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
+ parallel_process_termination();
+#endif
+ break;
+
+ case 10:
+ /* T='mantle_temp' for whole mantle + cold lithosphere TBL
+ + perturbations at some layers */
+
+ mantle_temperature = E->control.lith_age_mantle_temp;
+ constant_temperature_profile(E, mantle_temperature);
+ add_top_tbl(E, E->convection.half_space_age, mantle_temperature);
+ add_perturbations_at_all_layers(E);
+ break;
+
+ case 11:
+ /* T='mantle_temp' for whole mantle + hot CMB TBL
+ + perturbations at some layers */
+
+ mantle_temperature = E->control.lith_age_mantle_temp;
+ constant_temperature_profile(E, mantle_temperature);
+ add_bottom_tbl(E, E->convection.half_space_age, mantle_temperature);
+ add_perturbations_at_all_layers(E);
+ break;
+
+ case 12:
+ /* T='mantle_temp' for whole mantle + cold lithosphere TBL
+ + hot CMB TBL + perturbations at some layers */
+
+ mantle_temperature = E->control.lith_age_mantle_temp;
+ constant_temperature_profile(E, mantle_temperature);
+ add_top_tbl(E, E->convection.half_space_age, mantle_temperature);
+ add_bottom_tbl(E, E->convection.half_space_age, mantle_temperature);
+ add_perturbations_at_all_layers(E);
+ break;
+
+ case 100:
+ /* user-defined initial temperature goes here */
+ fprintf(stderr,"Need user definition for initial temperture: 'tic_method=%d'\n",
+ E->convection.tic_method);
+ parallel_process_termination();
+ break;
+
+ default:
+ /* unknown option */
+ fprintf(stderr,"Invalid value: 'tic_method=%d'\n", E->convection.tic_method);
+ parallel_process_termination();
+ break;
+ }
+
+ temperatures_conform_bcs(E);
+
+ /* debugging the code of expanding spherical harmonics */
+ /* debug_sphere_expansion(E);*/
+ return;
+}
+
+
Modified: mc/3D/CitcomS/trunk/lib/Regional_solver.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_solver.c 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/lib/Regional_solver.c 2008-11-07 23:32:20 UTC (rev 13273)
@@ -56,7 +56,6 @@
/* Version_dependent.c */
void regional_node_locations(struct All_variables *);
-void regional_construct_tic_from_input(struct All_variables *);
void regional_construct_boundary(struct All_variables *);
@@ -88,7 +87,6 @@
/* Version_dependent.c */
E->solver.node_locations = regional_node_locations;
- E->solver.construct_tic_from_input = regional_construct_tic_from_input;
E->solver.construct_boundary = regional_construct_boundary;
return;
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2008-11-07 23:32:20 UTC (rev 13273)
@@ -218,202 +218,6 @@
-void regional_construct_tic_from_input(struct All_variables *E)
-{
- double modified_plgndr_a(int, int, double);
- void temperatures_conform_bcs();
-
- int i, j ,k , kk, m, p, node, nodet;
- int nox, noy, noz, gnoz;
- double r1, f1, t1;
- int mm, ll;
- double con, temp;
-
- double theta_center;
- double fi_center;
- double r_center;
- double radius;
- double amp;
- double x_center,y_center,z_center;
- double theta,fi,r,x,y,z,distance;
-
- double tlen = M_PI / (E->control.theta_max - E->control.theta_min);
- double flen = M_PI / (E->control.fi_max - E->control.fi_min);
-
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- gnoz=E->mesh.noz;
-
- switch (E->convection.tic_method){
- case 0:
-
- /* 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;
-
- k = kk - E->lmesh.nzs + 1;
- if ( (k < 1) || (k >= noz) ) continue; /* if layer k is not inside this proc. */
- if (E->parallel.me_loc[1] == 0 && E->parallel.me_loc[2] == 0)
-
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=noy;i++)
- for(j=1;j<=nox;j++) {
- node=k+(j-1)*noz+(i-1)*nox*noz;
- t1 = (E->sx[m][1][node] - E->control.theta_min) * tlen;
- f1 = (E->sx[m][2][node] - E->control.fi_min) * flen;
-
- E->T[m][node] += con*cos(ll*t1)*cos(mm*f1);
-
- /*
- t1=E->sx[m][1][node];
- f1=E->sx[m][2][node];
- E->T[m][node] += con*modified_plgndr_a(ll,mm,t1)*cos(mm*f1);
- */
- }
- }
-
- break;
- case 1:
-
- /* set up a top thermal boundary layer */
- 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];
- temp = 0.2*(E->sphere.ro-r1) * 0.5/sqrt(E->convection.half_space_age/E->data.scalet);
- E->T[m][node] = E->control.TBCbotval*erf(temp);
- }
-
- break;
-
- case 2:
-
-
- theta_center = E->convection.blob_center[0];
- fi_center = E->convection.blob_center[1];
- r_center = E->convection.blob_center[2];
- radius = E->convection.blob_radius;
- amp = E->convection.blob_dT;
-
- fprintf(stderr,"center=%e %e %e radius=%e dT=%e\n",theta_center,fi_center,r_center,radius,amp);
- /* set up a thermal boundary layer 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];
- temp = 0.2*(E->sphere.ro-r1) * 0.5/sqrt(E->convection.half_space_age/E->data.scalet);
- E->T[m][node] = E->control.TBCbotval*erf(temp);
- }
-
- x_center = r_center * sin(fi_center) * cos(theta_center);
- y_center = r_center * sin(fi_center) * sin(theta_center);
- z_center = r_center * cos(fi_center);
-
- /* compute temperature field according to nodal coordinate */
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(k=1;k<=E->lmesh.noy;k++)
- for(j=1;j<=E->lmesh.nox;j++)
- for(i=1;i<=E->lmesh.noz;i++) {
- node = i + (j-1)*E->lmesh.noz
- + (k-1)*E->lmesh.noz*E->lmesh.nox;
-
- theta = E->sx[m][1][node];
- fi = E->sx[m][2][node];
- r = E->sx[m][3][node];
-
- distance = sqrt((theta - theta_center)*(theta - theta_center) +
- (fi - fi_center)*(fi - fi_center) +
- (r - r_center)*(r - r_center));
-
- if (distance < radius)
- E->T[m][node] += amp * exp(-1.0*distance/radius);
- }
- break;
- case 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));
- }
- }
- break;
- case 4: /* from grd files */
-#ifdef USE_GGRD
- ggrd_reg_temp_init(E);
-#else
- fprintf(stderr,"tic_method 4 only works for USE_GGRD compiled code\n");
- parallel_process_termination();
-#endif
- break;
-
- default: /* unknown option */
- fprintf(stderr,"Invalid value of 'tic_method'\n");
- parallel_process_termination();
- break;
- }
-
-
-
- temperatures_conform_bcs(E);
-
- return;
-}
-
-
/* setup boundary node and element arrays for bookkeeping */
void regional_construct_boundary( struct All_variables *E)
Modified: mc/3D/CitcomS/trunk/lib/solver.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/solver.h 2008-11-07 23:32:00 UTC (rev 13272)
+++ mc/3D/CitcomS/trunk/lib/solver.h 2008-11-07 23:32:20 UTC (rev 13273)
@@ -60,7 +60,6 @@
/* Version_dependent.c */
void (*node_locations)(struct All_variables *);
- void (*construct_tic_from_input)(struct All_variables *);
void (*construct_boundary)(struct All_variables *);
} solver;
More information about the CIG-COMMITS
mailing list