[cig-commits] r13270 - in mc/3D/CitcomS/trunk: CitcomS/Components/Sphere CitcomS/Solver examples/Full examples/Regional lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Nov 7 15:14:20 PST 2008
Author: tan2
Date: 2008-11-07 15:14:19 -0800 (Fri, 07 Nov 2008)
New Revision: 13270
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Sphere/Sphere.py
mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
mc/3D/CitcomS/trunk/examples/Full/input.sample
mc/3D/CitcomS/trunk/examples/Regional/input.sample
mc/3D/CitcomS/trunk/lib/Full_solver.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Instructions.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
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
* Added back parameters mgunitx/mgunity/mgunitz. This partly reverts r13256.
* Under pyre, uses mgunitx/mgunity/mgunitz and levels to compute nodex/nodey/nodez in multigrid solver.
* Merged regional_global_derived_values() and full_global_derived_values() to global_derived_values().
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Sphere/Sphere.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Sphere/Sphere.py 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Sphere/Sphere.py 2008-11-07 23:14:19 UTC (rev 13270)
@@ -84,6 +84,9 @@
nodex = pyre.inventory.int("nodex", default=9)
nodey = pyre.inventory.int("nodey", default=9)
nodez = pyre.inventory.int("nodez", default=9)
+ mgunitx = pyre.inventory.int("mgunitx", default=8)
+ mgunity = pyre.inventory.int("mgunity", default=8)
+ mgunitz = pyre.inventory.int("mgunitz", default=8)
levels = pyre.inventory.int("levels", default=1)
radius_outer = pyre.inventory.float("radius_outer", default=1.0)
Modified: mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py 2008-11-07 23:14:19 UTC (rev 13270)
@@ -251,9 +251,9 @@
Solver_set_properties(self.all_variables, self.inventory, stream)
inv = self.inventory
+ inv.vsolver.setProperties(stream)
inv.mesher.setProperties(stream)
inv.tsolver.setProperties(stream)
- inv.vsolver.setProperties(stream)
inv.bc.setProperties(stream)
inv.const.setProperties(stream)
Modified: mc/3D/CitcomS/trunk/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Full/input.sample 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/examples/Full/input.sample 2008-11-07 23:14:19 UTC (rev 13270)
@@ -32,6 +32,9 @@
nodex=9
nodey=9
nodez=9
+mgunitx=8
+mgunity=8
+mgunitz=8
levels=1
Modified: mc/3D/CitcomS/trunk/examples/Regional/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Regional/input.sample 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/examples/Regional/input.sample 2008-11-07 23:14:19 UTC (rev 13270)
@@ -33,6 +33,9 @@
nodex=17
nodey=17
nodez=9
+mgunitx=8
+mgunity=8
+mgunitz=8
levels=1
Modified: mc/3D/CitcomS/trunk/lib/Full_solver.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_solver.c 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/lib/Full_solver.c 2008-11-07 23:14:19 UTC (rev 13270)
@@ -55,7 +55,6 @@
void full_read_input_files_for_timesteps(struct All_variables *, int, int);
/* Version_dependent.c */
-void full_global_derived_values(struct All_variables *);
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 @@
E->solver.read_input_files_for_timesteps = full_read_input_files_for_timesteps;
/* Version_dependent.c */
- E->solver.global_derived_values = full_global_derived_values;
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;
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2008-11-07 23:14:19 UTC (rev 13270)
@@ -40,107 +40,7 @@
double theta_g(double , struct All_variables *);
#endif
-/* Setup global mesh parameters */
-void full_global_derived_values(E)
- struct All_variables *E;
-{
- int d,i,nox,noz,noy;
- E->mesh.levmax = E->mesh.levels-1;
- nox = E->mesh.mgunitx * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocx + 1;
- noy = E->mesh.mgunity * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocy + 1;
- noz = E->mesh.mgunitz * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocz + 1;
-
- if (E->control.NMULTIGRID||E->control.EMULTIGRID) {
- /* multigrid set up */
- E->mesh.levmax = E->mesh.levels-1;
- E->mesh.gridmax = E->mesh.levmax;
- E->mesh.nox = E->mesh.mgunitx * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocx + 1;
- E->mesh.noy = E->mesh.mgunity * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocy + 1;
- E->mesh.noz = E->mesh.mgunitz * (int) pow(2.0,((double)E->mesh.levmax))*E->parallel.nprocz + 1;
- } else {
- /* regular, conjugate gradient setup */
- if (nox!=E->mesh.nox || noy!=E->mesh.noy || noz!=E->mesh.noz) {
- if (E->parallel.me==0)
- fprintf(stderr,"inconsistent mesh for interpolation, quit the run\n");
- parallel_process_termination();
- }
- E->mesh.gridmax = E->mesh.levmax;
- E->mesh.gridmin = E->mesh.levmax;
- }
-
- if(E->mesh.nsd != 3)
- E->mesh.noy = 1;
-
- E->mesh.elx = E->mesh.nox-1;
- E->mesh.ely = E->mesh.noy-1;
- E->mesh.elz = E->mesh.noz-1;
-
- /* number of nodes, excluding overlaping nodes between processors */
- /* each cap has one row of nox and one row of noy overlapped, exclude these nodes.
- * nodes at north/south poles are exclued by all caps, include them by 2*noz*/
- E->mesh.nno = E->sphere.caps*(E->mesh.nox-1)*(E->mesh.noy-1)*E->mesh.noz + 2*E->mesh.noz;
-
- E->mesh.nel = E->sphere.caps*E->mesh.elx*E->mesh.elz*E->mesh.ely;
-
- E->mesh.nnov = E->mesh.nno;
-
- /* this is a rough estimate for global neq, a more accurate neq will
- be computed later. */
- E->mesh.neq = E->mesh.nnov*E->mesh.nsd;
-
- E->mesh.npno = E->mesh.nel;
- E->mesh.nsf = E->mesh.nox*E->mesh.noy;
-
- for(i=E->mesh.levmax;i>=E->mesh.levmin;i--) {
- if (E->control.NMULTIGRID||E->control.EMULTIGRID) {
- nox = E->mesh.mgunitx * (int) pow(2.0,(double)i)*E->parallel.nprocx + 1;
- noy = E->mesh.mgunity * (int) pow(2.0,(double)i)*E->parallel.nprocy + 1;
- noz = E->mesh.mgunitz * (int) pow(2.0,(double)i)*E->parallel.nprocz + 1;
- }
- else {
- noz = E->mesh.noz;
- nox = E->mesh.nox;
- noy = E->mesh.noy;
- /*if (i<E->mesh.levmax) noz=2;*/
- }
-
- E->mesh.ELX[i] = nox-1;
- E->mesh.ELY[i] = noy-1;
- E->mesh.ELZ[i] = noz-1;
- E->mesh.NNO[i] = E->sphere.caps * (nox-1) * (noy-1) * noz + 2 * noz;
- E->mesh.NEL[i] = E->sphere.caps * (nox-1) * (noz-1) * (noy-1);
- E->mesh.NPNO[i] = E->mesh.NEL[i] ;
- E->mesh.NOX[i] = nox;
- E->mesh.NOZ[i] = noz;
- E->mesh.NOY[i] = noy;
-
- E->mesh.NNOV[i] = E->mesh.NNO[i];
- E->mesh.NEQ[i] = E->mesh.nsd * E->mesh.NNOV[i] ;
- /* fprintf(stderr,"level=%d nox=%d noy=%d noz=%d %d %d %d %d %d %d %d %d %d %d %d\n",i,nox,noy,noz,E->mesh.ELX[i],E->mesh.ELY[i],E->mesh.ELZ[i],E->mesh.NNO[i],E->mesh.NEL[i],E->mesh.NPNO[i],E->mesh.NOX[i],E->mesh.NOZ[i],E->mesh.NOY[i],E->mesh.NNOV[i],E->mesh.NEQ[i]); */
- /* MPI_Barrier(E->parallel.world); */
- }
-
-
-
- /* Myr */
- E->data.scalet = (E->data.radius_km*1e3*E->data.radius_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
- /* cm/yr */
- E->data.scalev = (E->data.radius_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
- E->data.timedir = E->control.Atemp / fabs(E->control.Atemp);
-
- if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(stderr,"Problem has %d x %d x %d nodes per cap, %d nodes and %d elements in total\n",
- E->mesh.nox, E->mesh.noz, E->mesh.noy, E->mesh.nno, E->mesh.nel);
- fprintf(E->fp,"Problem has %d x %d x %d nodes per cap, %d nodes and %d elements in total\n",
- E->mesh.nox, E->mesh.noz, E->mesh.noy, E->mesh.nno, E->mesh.nel);
- }
-
- return;
-}
-
-
-
/* =================================================
rotate the mesh by a rotation matrix
=================================================*/
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-11-07 23:14:19 UTC (rev 13270)
@@ -84,6 +84,7 @@
void read_rayleigh_from_file(struct All_variables *);
void read_initial_settings(struct All_variables *);
void check_settings_consistency(struct All_variables *);
+void global_derived_values(struct All_variables *);
void initial_mesh_solver_setup(struct All_variables *E)
@@ -97,7 +98,7 @@
output_init(E);
(E->problem_derived_values)(E); /* call this before global_derived_ */
- (E->solver.global_derived_values)(E);
+ global_derived_values(E);
(E->solver.parallel_processor_setup)(E); /* get # of proc in x,y,z */
(E->solver.parallel_domain_decomp0)(E); /* get local nel, nno, elx, nox et al */
@@ -284,6 +285,7 @@
float tmp;
double ell_tmp;
int m=E->parallel.me;
+ double levmax;
/* first the problem type (defines subsequent behaviour) */
@@ -369,19 +371,28 @@
input_int("nprocy",&(E->parallel.nprocy),"1",m);
input_int("nprocz",&(E->parallel.nprocz),"1",m);
- input_int("nodex",&(E->mesh.nox),"essential",m);
- input_int("nodez",&(E->mesh.noz),"essential",m);
- input_int("nodey",&(E->mesh.noy),"essential",m);
+ if ( strcmp(E->control.SOLVER_TYPE,"cgrad") == 0) {
+ input_int("nodex",&(E->mesh.nox),"essential",m);
+ input_int("nodez",&(E->mesh.noz),"essential",m);
+ input_int("nodey",&(E->mesh.noy),"essential",m);
- input_int("levels",&(E->mesh.levels),"0",m);
+ E->mesh.mgunitx = E->mesh.nox - 1;
+ E->mesh.mgunity = E->mesh.noy - 1;
+ E->mesh.mgunitz = E->mesh.noz - 1;
+ E->mesh.levels = 1;
+ }
+ else {
+ input_int("mgunitx",&(E->mesh.mgunitx),"1",m);
+ input_int("mgunitz",&(E->mesh.mgunitz),"1",m);
+ input_int("mgunity",&(E->mesh.mgunity),"1",m);
- E->mesh.mgunitx = (E->mesh.nox - 1) / E->parallel.nprocx /
- (int) pow(2.0, E->mesh.levels - 1);
- E->mesh.mgunity = (E->mesh.noy - 1) / E->parallel.nprocy /
- (int) pow(2.0, E->mesh.levels - 1);
- E->mesh.mgunitz = (E->mesh.noz - 1) / E->parallel.nprocz /
- (int) pow(2.0, E->mesh.levels - 1);
+ input_int("levels",&(E->mesh.levels),"1",m);
+ levmax = E->mesh.levels - 1;
+ E->mesh.nox = E->mesh.mgunitx * (int) pow(2.0,levmax) * E->parallel.nprocx + 1;
+ E->mesh.noy = E->mesh.mgunity * (int) pow(2.0,levmax) * E->parallel.nprocy + 1;
+ E->mesh.noz = E->mesh.mgunitz * (int) pow(2.0,levmax) * E->parallel.nprocz + 1;
+ }
input_int("coor",&(E->control.coor),"0",m);
if(E->control.coor == 2){
@@ -754,6 +765,7 @@
/* Checking the consistency of input parameters */
void check_settings_consistency(struct All_variables *E)
{
+
if (strcmp(E->control.SOLVER_TYPE, "cgrad") == 0) {
/* conjugate gradient has only one level */
if(E->mesh.levels != 1)
@@ -771,6 +783,86 @@
}
+/* Setup global mesh parameters */
+void global_derived_values(struct All_variables *E)
+{
+ int d,i,nox,noz,noy;
+
+ E->mesh.levmax = E->mesh.levels-1;
+ E->mesh.gridmax = E->mesh.levmax;
+
+ E->mesh.elx = E->mesh.nox-1;
+ E->mesh.ely = E->mesh.noy-1;
+ E->mesh.elz = E->mesh.noz-1;
+
+ if(E->sphere.caps == 1) {
+ /* number of nodes, excluding overlaping nodes between processors */
+ E->mesh.nno = E->sphere.caps * E->mesh.nox * E->mesh.noy * E->mesh.noz;
+ }
+ else {
+ /* number of nodes, excluding overlaping nodes between processors */
+ /* each cap has one row of nox and one row of noy overlapped, exclude these nodes.
+ * nodes at north/south poles are exclued by all caps, include them by 2*noz*/
+ E->mesh.nno = E->sphere.caps * (E->mesh.nox-1) * (E->mesh.noy-1) * E->mesh.noz
+ + 2*E->mesh.noz;
+ }
+
+ E->mesh.nel = E->sphere.caps*E->mesh.elx*E->mesh.elz*E->mesh.ely;
+
+ E->mesh.nnov = E->mesh.nno;
+
+ /* this is a rough estimate for global neq, a more accurate neq will
+ be computed later. */
+ E->mesh.neq = E->mesh.nnov*E->mesh.nsd;
+
+ E->mesh.npno = E->mesh.nel;
+ E->mesh.nsf = E->mesh.nox*E->mesh.noy;
+
+ for(i=E->mesh.levmax;i>=E->mesh.levmin;i--) {
+ nox = E->mesh.mgunitx * (int) pow(2.0,(double)i)*E->parallel.nprocx + 1;
+ noy = E->mesh.mgunity * (int) pow(2.0,(double)i)*E->parallel.nprocy + 1;
+ noz = E->mesh.mgunitz * (int) pow(2.0,(double)i)*E->parallel.nprocz + 1;
+
+ E->mesh.ELX[i] = nox-1;
+ E->mesh.ELY[i] = noy-1;
+ E->mesh.ELZ[i] = noz-1;
+ if(E->sphere.caps == 1) {
+ E->mesh.NNO[i] = nox * noz * noy;
+ }
+ else {
+ E->mesh.NNO[i] = E->sphere.caps * (nox-1) * (noy-1) * noz + 2 * noz;
+ }
+ E->mesh.NEL[i] = E->sphere.caps * (nox-1) * (noz-1) * (noy-1);
+ E->mesh.NPNO[i] = E->mesh.NEL[i] ;
+ E->mesh.NOX[i] = nox;
+ E->mesh.NOZ[i] = noz;
+ E->mesh.NOY[i] = noy;
+
+ E->mesh.NNOV[i] = E->mesh.NNO[i];
+ E->mesh.NEQ[i] = E->mesh.nsd * E->mesh.NNOV[i] ;
+
+ }
+
+ /* Scaling from dimensionless units to Millions of years for input velocity
+ and time, timdir is the direction of time for advection. CPC 6/25/00 */
+
+ /* Myr */
+ E->data.scalet = (E->data.radius_km*1e3*E->data.radius_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
+ /* cm/yr */
+ E->data.scalev = (E->data.radius_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
+ E->data.timedir = E->control.Atemp / fabs(E->control.Atemp);
+
+
+ if(E->control.print_convergence && E->parallel.me==0) {
+ fprintf(stderr,"Problem has %d x %d x %d nodes per cap, %d nodes and %d elements in total\n",
+ E->mesh.nox, E->mesh.noz, E->mesh.noy, E->mesh.nno, E->mesh.nel);
+ fprintf(E->fp,"Problem has %d x %d x %d nodes per cap, %d nodes and %d elements in total\n",
+ E->mesh.nox, E->mesh.noz, E->mesh.noy, E->mesh.nno, E->mesh.nel);
+ }
+ return;
+}
+
+
/* ===================================
Functions which set up details
common to all problems follow ...
Modified: mc/3D/CitcomS/trunk/lib/Regional_solver.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_solver.c 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/lib/Regional_solver.c 2008-11-07 23:14:19 UTC (rev 13270)
@@ -55,7 +55,6 @@
void regional_read_input_files_for_timesteps(struct All_variables *, int, int);
/* Version_dependent.c */
-void regional_global_derived_values(struct All_variables *);
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 @@
E->solver.read_input_files_for_timesteps = regional_read_input_files_for_timesteps;
/* Version_dependent.c */
- E->solver.global_derived_values = regional_global_derived_values;
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;
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2008-11-07 23:14:19 UTC (rev 13270)
@@ -37,106 +37,7 @@
void ggrd_reg_temp_init(struct All_variables *);
#endif
-/* Setup global mesh parameters */
-void regional_global_derived_values(E)
- struct All_variables *E;
-{
- int d,i,nox,noz,noy;
- void parallel_process_termination();
-
-
- E->mesh.levmax = E->mesh.levels-1;
- nox = (E->mesh.mgunitx * (int) pow(2.0,((double)E->mesh.levmax)))*E->parallel.nprocx + 1;
- noy = (E->mesh.mgunity * (int) pow(2.0,((double)E->mesh.levmax)))*E->parallel.nprocy + 1;
-
- if (E->control.NMULTIGRID||E->control.EMULTIGRID) {
- E->mesh.levmax = E->mesh.levels-1;
- E->mesh.gridmax = E->mesh.levmax;
- E->mesh.nox = (E->mesh.mgunitx * (int) pow(2.0,((double)E->mesh.levmax)))*E->parallel.nprocx + 1;
- E->mesh.noy = (E->mesh.mgunity * (int) pow(2.0,((double)E->mesh.levmax)))*E->parallel.nprocy + 1;
- E->mesh.noz = (E->mesh.mgunitz * (int) pow(2.0,((double)E->mesh.levmax)))*E->parallel.nprocz + 1;
- }
- else {
- if (nox!=E->mesh.nox || noy!=E->mesh.noy) {
- if (E->parallel.me==0)
- fprintf(stderr,"inconsistent mesh for interpolation, quit the run\n");
- parallel_process_termination();
- }
- E->mesh.gridmax = E->mesh.levmax;
- E->mesh.gridmin = E->mesh.levmax;
- }
-
- if(E->mesh.nsd != 3)
- E->mesh.noy = 1;
-
- E->mesh.elx = E->mesh.nox-1;
- E->mesh.ely = E->mesh.noy-1;
- E->mesh.elz = E->mesh.noz-1;
-
- /* number of nodes, excluding overlaping nodes between processors */
- E->mesh.nno = E->sphere.caps*E->mesh.nox*E->mesh.noy*E->mesh.noz;
-
- E->mesh.nel = E->sphere.caps*E->mesh.elx*E->mesh.elz*E->mesh.ely;
-
- E->mesh.nnov = E->mesh.nno;
-
- /* this is a rough estimate for global neq, a more accurate neq will
- be computed later. */
- E->mesh.neq = E->mesh.nnov*E->mesh.nsd;
-
- E->mesh.npno = E->mesh.nel;
- E->mesh.nsf = E->mesh.nox*E->mesh.noy;
-
- for(i=E->mesh.levmax;i>=E->mesh.levmin;i--) {
- if (E->control.NMULTIGRID||E->control.EMULTIGRID)
- { nox = (E->mesh.mgunitx * (int) pow(2.0,(double)i))*E->parallel.nprocx + 1;
- noy = (E->mesh.mgunity * (int) pow(2.0,(double)i))*E->parallel.nprocy + 1;
- noz = (E->mesh.mgunitz * (int) pow(2.0,(double)i))*E->parallel.nprocz + 1;
- }
- else
- { noz = E->mesh.noz;
- nox = (E->mesh.mgunitx * (int) pow(2.0,(double)i))*E->parallel.nprocx + 1;
- noy = (E->mesh.mgunity * (int) pow(2.0,(double)i))*E->parallel.nprocy + 1;
- if (i<E->mesh.levmax) noz=2;
- }
-
- E->mesh.ELX[i] = nox-1;
- E->mesh.ELY[i] = noy-1;
- E->mesh.ELZ[i] = noz-1;
- E->mesh.NNO[i] = nox * noz * noy;
- E->mesh.NEL[i] = (nox-1) * (noz-1) * (noy-1);
- E->mesh.NPNO[i] = E->mesh.NEL[i] ;
- E->mesh.NOX[i] = nox;
- E->mesh.NOZ[i] = noz;
- E->mesh.NOY[i] = noy;
-
- E->mesh.NNOV[i] = E->mesh.NNO[i];
- E->mesh.NEQ[i] = E->mesh.nsd * E->mesh.NNOV[i] ;
-
- }
-
-/* Scaling from dimensionless units to Millions of years for input velocity
- and time, timdir is the direction of time for advection. CPC 6/25/00 */
-
- /* Myr */
- E->data.scalet = (E->data.radius_km*1e3*E->data.radius_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
- /* cm/yr */
- E->data.scalev = (E->data.radius_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
- E->data.timedir = E->control.Atemp / fabs(E->control.Atemp);
-
-
- if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(stderr,"Problem has %d x %d x %d nodes per cap, %d nodes and %d elements in total\n",
- E->mesh.nox, E->mesh.noz, E->mesh.noy, E->mesh.nno, E->mesh.nel);
- fprintf(E->fp,"Problem has %d x %d x %d nodes per cap, %d nodes and %d elements in total\n",
- E->mesh.nox, E->mesh.noz, E->mesh.noy, E->mesh.nno, E->mesh.nel);
- }
- return;
-}
-
-
-
/* =================================================
Standard node positions including mesh refinement
Modified: mc/3D/CitcomS/trunk/lib/solver.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/solver.h 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/lib/solver.h 2008-11-07 23:14:19 UTC (rev 13270)
@@ -59,7 +59,6 @@
void (*read_input_files_for_timesteps)(struct All_variables *, int, int);
/* Version_dependent.c */
- void (*global_derived_values)(struct All_variables *);
void (*node_locations)(struct All_variables *);
void (*construct_tic_from_input)(struct All_variables *);
void (*construct_boundary)(struct All_variables *);
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2008-11-07 05:47:46 UTC (rev 13269)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2008-11-07 23:14:19 UTC (rev 13270)
@@ -545,18 +545,30 @@
getFloatVectorProperty(properties, "coor_refine", E->control.coor_refine, 4, fp);
getStringProperty(properties, "coor_file", E->control.coor_file, fp);
- getIntProperty(properties, "nodex", E->mesh.nox, fp);
- getIntProperty(properties, "nodey", E->mesh.noy, fp);
- getIntProperty(properties, "nodez", E->mesh.noz, fp);
- getIntProperty(properties, "levels", E->mesh.levels, fp);
+ if ( strcmp(E->control.SOLVER_TYPE,"cgrad") == 0) {
+ getIntProperty(properties, "nodex", E->mesh.nox, fp);
+ getIntProperty(properties, "nodey", E->mesh.noy, fp);
+ getIntProperty(properties, "nodez", E->mesh.noz, fp);
- E->mesh.mgunitx = (E->mesh.nox - 1) / E->parallel.nprocx /
- (int) pow(2.0, E->mesh.levels - 1);
- E->mesh.mgunity = (E->mesh.noy - 1) / E->parallel.nprocy /
- (int) pow(2.0, E->mesh.levels - 1);
- E->mesh.mgunitz = (E->mesh.noz - 1) / E->parallel.nprocz /
- (int) pow(2.0, E->mesh.levels - 1);
+ E->mesh.mgunitx = E->mesh.nox - 1;
+ E->mesh.mgunity = E->mesh.noy - 1;
+ E->mesh.mgunitz = E->mesh.noz - 1;
+ E->mesh.levels = 1;
+ }
+ else {
+ double levmax;
+ getIntProperty(properties, "mgunitx", E->mesh.mgunitx, fp);
+ getIntProperty(properties, "mgunity", E->mesh.mgunity, fp);
+ getIntProperty(properties, "mgunitz", E->mesh.mgunitz, fp);
+ getIntProperty(properties, "levels", E->mesh.levels, fp);
+
+ levmax = E->mesh.levels - 1;
+ E->mesh.nox = E->mesh.mgunitx * (int) pow(2.0,levmax) * E->parallel.nprocx + 1;
+ E->mesh.noy = E->mesh.mgunity * (int) pow(2.0,levmax) * E->parallel.nprocy + 1;
+ E->mesh.noz = E->mesh.mgunitz * (int) pow(2.0,levmax) * E->parallel.nprocz + 1;
+ }
+
if (E->parallel.nprocxy == 12) {
if (E->mesh.nox != E->mesh.noy) {
char errmsg[] = "!!!! nodex must equal to nodey";
More information about the CIG-COMMITS
mailing list