[cig-commits] r14786 - in mc/3D/CitcomS/branches/cxx: bin lib
leif at geodynamics.org
leif at geodynamics.org
Wed Apr 22 12:42:54 PDT 2009
Author: leif
Date: 2009-04-22 12:42:53 -0700 (Wed, 22 Apr 2009)
New Revision: 14786
Modified:
mc/3D/CitcomS/branches/cxx/bin/Citcom.cc
mc/3D/CitcomS/branches/cxx/lib/Checkpoints.cc
mc/3D/CitcomS/branches/cxx/lib/Composition_related.cc
mc/3D/CitcomS/branches/cxx/lib/Determine_net_rotation.cc
mc/3D/CitcomS/branches/cxx/lib/Full_tracer_advection.cc
mc/3D/CitcomS/branches/cxx/lib/Instructions.cc
mc/3D/CitcomS/branches/cxx/lib/Regional_tracer_advection.cc
mc/3D/CitcomS/branches/cxx/lib/Tracer_setup.cc
Log:
Merged r14108:14275 from trunk.
Modified: mc/3D/CitcomS/branches/cxx/bin/Citcom.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/bin/Citcom.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/bin/Citcom.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -37,6 +37,7 @@
#include "parallel_util.h"
#include "checkpoints.h"
+#include "composition_related.h"
#include "drive_solvers.h"
#include "instructions.h"
#include "interuption.h"
@@ -52,7 +53,7 @@
{ /* Functions called by main*/
float cpu_time_on_vp_it;
- int cpu_total_seconds,k, *temp;
+ int cpu_total_seconds,k,need_init_sol;
double time,initial_time,start_time;
struct All_variables *E;
@@ -104,12 +105,28 @@
if (E->control.restart) {
/* the initial condition is from previous checkpoint */
read_checkpoint(E);
+ if(E->control.tracer && (E->trace.ic_method_for_flavors == 99)){
+ /*
+ if ggrd tracer input is selected, this will override
+ existing tracers, or allow addition of tracers after
+ restart from a thermal model
+
+ */
+ initialize_tracers(E);
+ if (E->composition.on)
+ init_composition(E);
+ need_init_sol = 1;
+ }else
+ need_init_sol = 0;
}
else {
/* regular init, or read T from file only */
initial_conditions(E);
-
+ need_init_sol = 1; /* */
+ }
+ if(need_init_sol){
+ /* find first solution */
if(E->control.pseudo_free_surf) {
if(E->mesh.topvbc == 2)
general_stokes_solver_pseudo_surf(E);
Modified: mc/3D/CitcomS/branches/cxx/lib/Checkpoints.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Checkpoints.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Checkpoints.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -124,10 +124,14 @@
/* read tracer/composition information in the checkpoint file */
if(E->control.tracer) {
+ if(E->trace.ic_method_for_flavors == 99){
+ if(E->parallel.me == 0)
+ fprintf(stderr,"ic_method_for_flavors = 99 will override checkpoint restart\n");
+ }else{
read_tracer_checkpoint(E, fp);
-
if(E->composition.on)
read_composition_checkpoint(E, fp);
+ }
}
fclose(fp);
Modified: mc/3D/CitcomS/branches/cxx/lib/Composition_related.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Composition_related.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Composition_related.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -84,6 +84,7 @@
/* compositional rheology */
+ /* what was this about? there is a CDEPV for compositional rheology TWB */
/* icompositional_rheology=0 (off) */
/* icompositional_rheology=1 (on) */
Modified: mc/3D/CitcomS/branches/cxx/lib/Determine_net_rotation.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Determine_net_rotation.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Determine_net_rotation.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -53,8 +53,10 @@
static double determine_netr_tp(float r, float theta, float phi, float velt, float velp, int mode, double *c9, double *omega);
-static void hc_ludcmp_3x3(double a[3][3], int *indx);
-static void hc_lubksb_3x3(double a[3][3], int *indx, double *b);
+#ifdef USE_GGRD
+static void hc_ludcmp_3x3(HC_PREC [3][3],int,int *);
+static void hc_lubksb_3x3(HC_PREC [3][3],int,int *,HC_PREC *);
+#endif
/*
@@ -221,6 +223,7 @@
{
float coslat,coslon,sinlat,sinlon,rx,ry,rz,rate,rzu,a,b,c,d,e,f;
int i,j,ind[3];
+ static int n3 = 3;
double amp,coef[3][3];
switch(mode){
case 0: /* initialize */
@@ -277,9 +280,11 @@
omega[1] = c9[7];
omega[2] = c9[8];
- /* solve solution*/
- hc_ludcmp_3x3(coef,ind);
- hc_lubksb_3x3(coef,ind,omega);
+ /* solve */
+#ifdef USE_GGRD
+ hc_ludcmp_3x3(coef,n3,ind);
+ hc_lubksb_3x3(coef,n3,ind,omega);
+#endif
amp = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
break;
default:
@@ -370,45 +375,46 @@
// AND MAILING.
//
+#ifdef USE_GGRD
+#define NR_TINY 1.0e-20;
+/*
+ matrix is always 3 x 3 , solution is for full system for n == 3,
+ for upper 2 x 2 only for n = 2
-/*
-
-matrix solvers from numerical recipes
-
*/
-#define NR_TINY 1.0e-20;
-
-static void hc_ludcmp_3x3(double a[3][3],int *indx)
+static void hc_ludcmp_3x3(HC_PREC a[3][3],int n,int *indx)
{
int i,imax=0,j,k;
- double big,dum,sum,temp;
- double vv[3];
-
- for (i=0;i < 3;i++) {
+ HC_PREC big,dum,sum,temp;
+ HC_PREC vv[3];
+
+ for (i=0;i < n;i++) {
big=0.0;
- for (j=0;j < 3;j++)
- if ((temp = fabs(a[i][j])) > big)
+ for (j=0;j < n;j++)
+ if ((temp = fabs(a[i][j])) > big)
big=temp;
- if (fabs(big) < 5e-15) {
+ if (fabs(big) < HC_EPS_PREC) {
fprintf(stderr,"hc_ludcmp_3x3: singular matrix in routine, big: %g\n",
big);
- //hc_print_3x3(a,stderr);
- for(j=0;j<3;j++)
- fprintf(stderr,"%g %g %g\n",a[j][0],a[j][1],a[j][2]);
- parallel_process_termination();
+ for(j=0;j <n;j++){
+ for(k=0;k<n;k++)
+ fprintf(stderr,"%g ",a[j][k]);
+ fprintf(stderr,"\n");
+ }
+ exit(-1);
}
vv[i]=1.0/big;
}
- for (j=0;j < 3;j++) {
+ for (j=0;j < n;j++) {
for (i=0;i < j;i++) {
sum = a[i][j];
- for (k=0;k < i;k++)
+ for (k=0;k < i;k++)
sum -= a[i][k] * a[k][j];
a[i][j]=sum;
}
big=0.0;
- for (i=j;i < 3;i++) {
+ for (i=j;i < n;i++) {
sum=a[i][j];
for (k=0;k < j;k++)
sum -= a[i][k] * a[k][j];
@@ -419,7 +425,7 @@
}
}
if (j != imax) {
- for (k=0;k < 3;k++) {
+ for (k=0;k < n;k++) {
dum = a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
@@ -427,36 +433,41 @@
vv[imax]=vv[j];
}
indx[j]=imax;
- if (fabs(a[j][j]) < 5e-15)
+ if (fabs(a[j][j]) < HC_EPS_PREC)
a[j][j] = NR_TINY;
if (j != 2) {
dum=1.0/(a[j][j]);
- for (i=j+1;i < 3;i++)
+ for (i=j+1;i < n;i++)
a[i][j] *= dum;
}
}
}
+
#undef NR_TINY
-static void hc_lubksb_3x3(double a[3][3], int *indx, double *b)
+static void hc_lubksb_3x3(HC_PREC a[3][3], int n,int *indx, HC_PREC *b)
{
int i,ii=0,ip,j;
- double sum;
- for (i=0;i < 3;i++) {
+ HC_PREC sum;
+ int nm1;
+ nm1 = n - 1;
+ for (i=0;i < n;i++) {
ip = indx[i];
sum = b[ip];
b[ip]=b[i];
if (ii)
- for (j=ii-1;j <= i-1;j++)
+ for (j=ii-1;j <= i-1;j++)
sum -= a[i][j]*b[j];
- else if (fabs(sum) > 5e-15)
+ else if (fabs(sum)>HC_EPS_PREC)
ii = i+1;
b[i]=sum;
}
- for (i=2;i>=0;i--) {
+ for (i=nm1;i>=0;i--) {
sum=b[i];
- for (j=i+1;j < 3;j++)
+ for (j=i+1;j < n;j++)
sum -= a[i][j]*b[j];
b[i] = sum/a[i][i];
}
}
+#endif
+
Modified: mc/3D/CitcomS/branches/cxx/lib/Full_tracer_advection.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Full_tracer_advection.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Full_tracer_advection.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -1920,12 +1920,14 @@
if (E->trace.nflavors && E->trace.ic_method==0) {
fprintf(E->trace.fpt,"Initialized tracer flavors by: %d\n", E->trace.ic_method_for_flavors);
if (E->trace.ic_method_for_flavors == 0) {
+ /* default mode 0 */
fprintf(E->trace.fpt,"Layered tracer flavors\n");
for (i=0; i<E->trace.nflavors-1; i++)
fprintf(E->trace.fpt,"Interface Height: %d %f\n",i,E->trace.z_interface[i]);
}
#ifdef USE_GGRD
- else if(E->trace.ic_method_for_flavors == 1) {
+ else if((E->trace.ic_method_for_flavors == 1)||(E->trace.ic_method_for_flavors == 99)) {
+ /* ggrd modes 1 and 99 (99 is override for restart) */
fprintf(E->trace.fpt,"netcdf grd assigned tracer flavors\n");
fprintf(E->trace.fpt,"file: %s top %i layeres\n",E->trace.ggrd_file,
E->trace.ggrd_layers);
Modified: mc/3D/CitcomS/branches/cxx/lib/Instructions.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Instructions.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Instructions.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -175,6 +175,7 @@
if(E->control.tracer) {
tracer_initial_settings(E);
(E->problem_tracer_setup)(E);
+ if(chatty)fprintf(stderr,"tracer setup done\n");
}
if(chatty)fprintf(stderr,"initial_mesh_solver_setup done\n");
}
Modified: mc/3D/CitcomS/branches/cxx/lib/Regional_tracer_advection.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Regional_tracer_advection.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Regional_tracer_advection.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -182,6 +182,16 @@
for (i=0; i<E->trace.nflavors-1; i++)
fprintf(E->trace.fpt,"Interface Height: %d %f\n",i,E->trace.z_interface[i]);
}
+#ifdef USE_GGRD
+ else if((E->trace.ic_method_for_flavors == 1)||(E->trace.ic_method_for_flavors == 99)) {
+ /* ggrd modes 1 and 99 (99 is override for restart) */
+ fprintf(stderr,"ggrd regional flavors not implemented\n");
+ fprintf(E->trace.fpt,"ggrd not implemented et for regional, flavor method= %d\n",
+ E->trace.ic_method_for_flavors);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+#endif
else {
fprintf(E->trace.fpt,"Sorry-This IC methods for Flavors are Unavailable %d\n",E->trace.ic_method_for_flavors);
fflush(E->trace.fpt);
Modified: mc/3D/CitcomS/branches/cxx/lib/Tracer_setup.cc
===================================================================
--- mc/3D/CitcomS/branches/cxx/lib/Tracer_setup.cc 2009-04-22 00:59:20 UTC (rev 14785)
+++ mc/3D/CitcomS/branches/cxx/lib/Tracer_setup.cc 2009-04-22 19:42:53 UTC (rev 14786)
@@ -136,6 +136,13 @@
* either using absolute method or ratio method. */
input_int("tracer_flavors",&(E->trace.nflavors),"0,0,nomax",m);
+ /* 0: default from layers
+ 1: from netcdf grds
+
+
+ 99: from grds, overriding checkpoints during restart
+ (1 and 99 require ggrd)
+ */
input_int("ic_method_for_flavors",
&(E->trace.ic_method_for_flavors),"0,0,nomax",m);
@@ -143,7 +150,9 @@
if (E->trace.nflavors > 1) {
switch(E->trace.ic_method_for_flavors){
- case 0: /* layer */
+ /* default method */
+ case 0:
+ /* flavors initialized from layers */
E->trace.z_interface = (double*) malloc((E->trace.nflavors-1)
*sizeof(double));
for(i=0; i<E->trace.nflavors-1; i++)
@@ -152,12 +161,21 @@
input_double_vector("z_interface", E->trace.nflavors-1,
E->trace.z_interface, m);
break;
- case 1: /* from grid in top n materials */
- input_string("ictracer_grd_file",E->trace.ggrd_file,"",m); /* file from which to read */
- input_int("ictracer_grd_layers",&(E->trace.ggrd_layers),"2",m); /* which top layers to use */
- break;
+ /*
+ two grd init method, second will override restart
+ */
+#ifdef USE_GGRD
+ case 1:
+ case 99: /* will override restart */
+ /* from grid in top n materials, this will override
+ the checkpoint input */
+ input_string("ictracer_grd_file",E->trace.ggrd_file,"",m); /* file from which to read */
+ input_int("ictracer_grd_layers",&(E->trace.ggrd_layers),"2",m); /* which top layers to use */
+ break;
+
+#endif
default:
- fprintf(stderr,"ic_method_for_flavors %i undefined\n",E->trace.ic_method_for_flavors);
+ fprintf(stderr,"ic_method_for_flavors %i undefined (1 and 99 only for ggrd mode)\n",E->trace.ic_method_for_flavors);
parallel_process_termination();
break;
}
@@ -1207,6 +1225,7 @@
break;
case 1: /* from grd in top n layers */
+ case 99: /* (will override restart) */
#ifndef USE_GGRD
fprintf(stderr,"ic_method_for_flavors %i requires the ggrd routines from hc, -DUSE_GGRD\n",
E->trace.ic_method_for_flavors);
More information about the CIG-COMMITS
mailing list