[cig-commits] r18591 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Sat Jun 11 19:11:10 PDT 2011
Author: becker
Date: 2011-06-11 19:11:10 -0700 (Sat, 11 Jun 2011)
New Revision: 18591
Modified:
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
Log:
Adjusted ggrd temperature and composition input such that different settings can be used for slices.
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2011-06-11 00:25:45 UTC (rev 18590)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2011-06-12 02:11:10 UTC (rev 18591)
@@ -206,7 +206,7 @@
report(E, "convection, initial temperature");
#ifdef USE_GGRD
- convection_initial_temperature_ggrd(E);
+ convection_initial_temperature_and_comp_ggrd(E);
#else
convection_initial_temperature(E);
#endif
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2011-06-11 00:25:45 UTC (rev 18590)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2011-06-12 02:11:10 UTC (rev 18591)
@@ -49,7 +49,7 @@
Initialization of fields .....
=============================== */
-void convection_initial_temperature_ggrd(struct All_variables *E)
+void convection_initial_temperature_and_comp_ggrd(struct All_variables *E)
{
int ll, mm, i, j, k, p, node, ii,slice,hit;
double temp, temp1, temp2, temp3, base, radius, radius2;
@@ -90,216 +90,218 @@
top_t = (E->mesh.toptbc) ? E->control.TBCtopval : 0.0;
for(i=1;i<=E->lmesh.nno;i++) {
- /* init as zeros */
+ /*
+ init both composition and temperature as zeros
+ */
E->T[i] = E->C[i] = 0.0;
}
- if(!E->control.restart)
- { /*
-
- regular init
-
- */
-
- if(E->control.ggrd.use_temp){ /* read T init from grid */
- if(E->parallel.me==0)
- fprintf(stderr,"convection_initial_temperature: using GMT grd files for temperatures\n");
- /*
+ if(!E->control.restart){
+ /*
+
+ regular init
+
+ */
+
+ if(E->control.ggrd.use_temp){ /* read T init from grid */
+ if(E->parallel.me==0)
+ fprintf(stderr,"convection_initial_temperature: using GMT grd files for temperatures\n");
+ /*
- read in tempeatures/density from GMT grd files
+ read in tempeatures/density from GMT grd files
- */
- if(E->parallel.me > 0){
- /*
- wait for the previous processor
- */
- mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
- MPI_COMM_WORLD, &mpi_stat);
- }
- if(E->parallel.me == 0)
- fprintf(stderr,"Init: initializing PREM and ggrd files\n");
- if(E->control.ggrd.temp.scale_with_prem){
- if(!E->control.Rsphere)
- myerror("cannot use PREM with Cartesian setting",E);
- /* initialize PREM */
- if(prem_read_model(PREM_MODEL_FILE,&E->control.ggrd.temp.prem,
- (E->parallel.me==0)))
- myerror("prem initialization",E);
- }
+ */
+ if(E->parallel.me > 0){
/*
- initialize the GMT grid files
+ wait for the previous processor
*/
- if(E->control.ggrd_slab_slice){ /* only a few slices of x - depth */
- if(E->control.ggrd_slab_slice == 1){
- if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.temp.gfile,
- char_dummy,"",E->control.ggrd_ss_grd,
+ mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
+ MPI_COMM_WORLD, &mpi_stat);
+ }
+ if(E->parallel.me == 0)
+ fprintf(stderr,"Init: initializing PREM and ggrd files\n");
+ if(E->control.ggrd.temp.scale_with_prem){
+ if(!E->control.Rsphere)
+ myerror("cannot use PREM with Cartesian setting",E);
+ /* initialize PREM */
+ if(prem_read_model(PREM_MODEL_FILE,&E->control.ggrd.temp.prem,
+ (E->parallel.me==0)))
+ myerror("prem initialization",E);
+ }
+ /*
+ initialize the GMT grid files for temperatures
+ */
+ if(E->control.ggrd_t_slab_slice){ /* only a few slices of x - depth */
+ if(E->control.ggrd_t_slab_slice == 1){
+ if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd.temp.gfile,
+ char_dummy,"",E->control.ggrd_ss_grd,
+ (E->parallel.me==0),
+ FALSE,FALSE))
+ myerror("grdtrack x-z init error",E);
+ }else{ /* read several slab slices */
+ for(slice=0;slice < E->control.ggrd_t_slab_slice;slice++){
+ sprintf(pfile,"%s.%i.grd",E->control.ggrd.temp.gfile,slice+1);
+ if(ggrd_grdtrack_init_general(FALSE,pfile,
+ char_dummy,"",(E->control.ggrd_ss_grd+slice),
(E->parallel.me==0),
FALSE,FALSE))
myerror("grdtrack x-z init error",E);
- }else{ /* several slab slices */
- for(slice=0;slice < E->control.ggrd_slab_slice;slice++){
- sprintf(pfile,"%s.%i.grd",E->control.ggrd.temp.gfile,slice+1);
- if(ggrd_grdtrack_init_general(FALSE,pfile,
- char_dummy,"",(E->control.ggrd_ss_grd+slice),
- (E->parallel.me==0),
- FALSE,FALSE))
- myerror("grdtrack x-z init error",E);
- }
}
- }else{ /* 3-D */
- if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.temp.gfile,
- E->control.ggrd.temp.dfile,"",
- E->control.ggrd.temp.d,(E->parallel.me==0),
- FALSE,FALSE))
- myerror("grdtrack 3-D init error",E);
}
+ }else{ /* 3-D */
+ if(ggrd_grdtrack_init_general(TRUE,E->control.ggrd.temp.gfile,
+ E->control.ggrd.temp.dfile,"",
+ E->control.ggrd.temp.d,(E->parallel.me==0),
+ FALSE,FALSE))
+ myerror("grdtrack 3-D init error",E);
+ }
- if(E->parallel.me < E->parallel.nproc-1){
- /*
- tell the next processor to go ahead
- */
- mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
- (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
- }else{
- fprintf(stderr,"convection_initial_temperature: last processor done with temp grd init\n");
- }
+ if(E->parallel.me < E->parallel.nproc-1){
/*
+ tell the next processor to go ahead
+ */
+ mpi_rc = MPI_Send(&mpi_success_message, 1, MPI_INT,
+ (E->parallel.me+1), mpi_tag, MPI_COMM_WORLD);
+ }else{
+ fprintf(stderr,"convection_initial_temperature: last processor done with temp grd init\n");
+ }
+ /*
- end MPI synchornization part
+ end MPI synchornization part
- */
- /*
+ */
+ /*
- interpolate densities to temperature given PREM variations
+ interpolate densities to temperature given PREM variations
- */
+ */
- tmean = (top_t + bot_t)/2.0 + E->control.ggrd.temp.offset;
- if(E->parallel.me == 0)
- fprintf(stderr,"tinit: top: %g bot: %g mean: %g scale: %g PREM %i\n",
- top_t,bot_t,tmean,E->control.ggrd.temp.scale,
- E->control.ggrd.temp.scale_with_prem);
- 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; /* offset */
- if(E->control.ggrd_slab_slice){
- /*
+ tmean = (top_t + bot_t)/2.0 + E->control.ggrd.temp.offset;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"tinit: top: %g bot: %g mean: %g scale: %g PREM %i\n",
+ top_t,bot_t,tmean,E->control.ggrd.temp.scale,
+ E->control.ggrd.temp.scale_with_prem);
+ 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; /* offset */
+ if(E->control.ggrd_t_slab_slice){
+ /*
- slab slice input
+ slab slice input
- */
- for(slice=hit=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
- if(E->control.Rsphere) {
- if(in_slab_slice(E->SX[1][node],slice,E)){
- /* spherical interpolation for a slice ! */
- ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
- (double)E->SX[1][node],
- (E->control.ggrd_ss_grd+slice),&tadd,
- FALSE);
- hit=1;
- }
- }else{ /* cartesian interpolation */
- if(in_slab_slice(E->X[2][node],slice,E)){
- ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
- (double)E->X[3][node],
- (E->control.ggrd_ss_grd+slice),&tadd,
- FALSE);
- hit = 1;
- }
+ */
+ for(slice=hit=0;(!hit) && (slice < E->control.ggrd_t_slab_slice);slice++){
+ if(E->control.Rsphere) {
+ if(in_slab_slice(E->SX[1][node],slice,E, E->control.ggrd_t_slab_theta_bound,E->control.ggrd_t_slab_slice)){
+ /* spherical interpolation for a slice ! */
+ ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
+ (double)E->SX[1][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit=1;
}
+ }else{ /* cartesian interpolation */
+ if(in_slab_slice(E->X[2][node],slice,E, E->control.ggrd_t_slab_theta_bound,E->control.ggrd_t_slab_slice)){
+ ggrd_grdtrack_interpolate_xy((double)E->X[1][node],
+ (double)E->X[3][node],
+ (E->control.ggrd_ss_grd+slice),&tadd,
+ FALSE);
+ hit = 1;
+ }
}
- if(!hit)
- if(((E->control.Rsphere) && (E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))||
- (E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))
- tadd = E->control.TBCtopval;
- else
- tadd = 1.0;
- /* end slice part */
- }else{
- /*
+ }
+ if(!hit)
+ if(((E->control.Rsphere) && (E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))||
+ (E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))
+ tadd = E->control.TBCtopval;
+ else
+ tadd = 1.0;
+ /* end slice part */
+ }else{
+ /*
- 3-D grid model input
+ 3-D grid model input
- */
- if(E->control.Rsphere) /* spherical interpolation */
- ggrd_grdtrack_interpolate_rtp((double)E->SX[3][node],
- (double)E->SX[1][node],
- (double)E->SX[2][node],
- E->control.ggrd.temp.d,&tadd,
- FALSE,FALSE);
- else{ /* cartesian interpolation */
- ggrd_grdtrack_interpolate_xyz((double)E->X[1][node],
- (double)E->X[2][node],
- (double)E->X[3][node],
- E->control.ggrd.temp.d,&tadd,
- FALSE);
- }
+ */
+ if(E->control.Rsphere) /* spherical interpolation */
+ ggrd_grdtrack_interpolate_rtp((double)E->SX[3][node],
+ (double)E->SX[1][node],
+ (double)E->SX[2][node],
+ E->control.ggrd.temp.d,&tadd,
+ FALSE,FALSE);
+ else{ /* cartesian interpolation */
+ ggrd_grdtrack_interpolate_xyz((double)E->X[1][node],
+ (double)E->X[2][node],
+ (double)E->X[3][node],
+ E->control.ggrd.temp.d,&tadd,
+ FALSE);
}
- if(E->control.ggrd.temp.scale_with_prem){ /* only works for spherical! */
- /*
- get the PREM density at r for additional scaling
- */
- prem_get_rho(&rho_prem,(double)E->SX[3][node],&E->control.ggrd.temp.prem);
- if(rho_prem < 3200.0)
- rho_prem = 3200.0; /* we don't want the viscosity of water */
- /*
- assign temperature
- */
- E->T[node] = tmean + tadd * E->control.ggrd.temp.scale *
- rho_prem / E->data.density;
- }else{
- /* no PREM scaling */
- E->T[node] = tmean + tadd * E->control.ggrd.temp.scale;
-
- //fprintf(stderr,"z: %11g mean: %11g tadd: %11g scale: %11g T: %11g\n", E->X[3][node],tmean, tadd ,E->control.ggrd_tinit_scale, E->T[node] );
-
- }
+ }
+ if(E->control.ggrd.temp.scale_with_prem){ /* only works for spherical! */
/*
- if we're on the surface or bottom, reassign T and
- temperature boundary condition if toptbc or bottbsc == 2
+ get the PREM density at r for additional scaling
*/
- if(E->control.Rsphere){
- if((E->mesh.toptbc == 2) && (E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1])){
- E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
- }
- if((E->mesh.bottbc == 2) && (E->SX[3][node] == E->segment.zzlayer[0])){
- E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
- }
- }else{
- if((E->mesh.toptbc == 2) && (E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))
- E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
- if((E->mesh.bottbc == 2) && (E->X[3][node] == E->segment.zzlayer[0]))
- E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
- }
+ prem_get_rho(&rho_prem,(double)E->SX[3][node],&E->control.ggrd.temp.prem);
+ if(rho_prem < 3200.0)
+ rho_prem = 3200.0; /* we don't want the viscosity of water */
/*
- boundary condition flags
+ assign temperature
*/
- if(!setflag)
- E->node[node] = E->node[node] | (INTX | INTZ | INTY);
- } /* end node loop */
- setflag=1;
- /* free the structure, not needed anymore */
- ggrd_grdtrack_free_gstruc(E->control.ggrd.temp.d);
- /*
- end temperature from GMT grd init
- */
- /* end grid init branch */
- } else {
+ E->T[node] = tmean + tadd * E->control.ggrd.temp.scale *
+ rho_prem / E->data.density;
+ }else{
+ /* no PREM scaling */
+ E->T[node] = tmean + tadd * E->control.ggrd.temp.scale;
+
+ //fprintf(stderr,"z: %11g mean: %11g tadd: %11g scale: %11g T: %11g\n", E->X[3][node],tmean, tadd ,E->control.ggrd_tinit_scale, E->T[node] );
+ }
+ /*
+ if we're on the surface or bottom, reassign T and
+ temperature boundary condition if toptbc or bottbsc == 2
+ */
+ if(E->control.Rsphere){
+ if((E->mesh.toptbc == 2) && (E->SX[3][node] == E->segment.zzlayer[E->segment.zlayers-1])){
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ }
+ if((E->mesh.bottbc == 2) && (E->SX[3][node] == E->segment.zzlayer[0])){
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ }
+ }else{
+ if((E->mesh.toptbc == 2) && (E->X[3][node] == E->segment.zzlayer[E->segment.zlayers-1]))
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ if((E->mesh.bottbc == 2) && (E->X[3][node] == E->segment.zzlayer[0]))
+ E->TB[1][node] = E->TB[2][node] = E->TB[3][node] = E->T[node];
+ }
+ /*
+ boundary condition flags
+ */
+ if(!setflag)
+ E->node[node] = E->node[node] | (INTX | INTZ | INTY);
+ } /* end node loop */
+ setflag=1;
+ /* free the structure, not needed anymore */
+ ggrd_grdtrack_free_gstruc(E->control.ggrd.temp.d);
+ /*
+ end temperature from GMT grd init
+ */
+ /* end grid init branch */
+ } else {
- /*
+
+ /*
- add a linear profile and potentially perturbations to the temperature
+ add a linear profile and potentially perturbations to the temperature
- */
+ */
- if(E->control.CART3D)
- {
+ if(E->control.CART3D)
+ {
for(i = 1; i <= noy; i++)
for(j = 1; j <= nox; j++)
@@ -367,45 +369,45 @@
setflag =1;
} // end for if else if of geometry
- } /* end perturnbation branch */
- /*
+ } /* end perturnbation branch */
+ /*
- now deal with composition
+ now deal with composition
- */
- if(E->control.ggrd.use_comp){ /* read composition init from grid */
- if(!E->control.composition)
- myerror("comp grd init but no composition control set!?",E);
+ */
+ if(E->control.ggrd.use_comp){ /* read composition init from grid */
+ if(!E->control.composition)
+ myerror("comp grd init but no composition control set!?",E);
- ggrd_deal_with_composition_input(E, 1); /* compositional (density) */
+ ggrd_deal_with_composition_input(E, 1); /* compositional (density) */
- /* flavors are dealt with based on markers, and then assigned
- to nodes */
+ /* flavors are dealt with based on markers, and then assigned
+ to nodes */
- //if(E->tracers_add_flavors) /* flavors */
- //ggrd_deal_with_composition_input(E, 0);
+ //if(E->tracers_add_flavors) /* flavors */
+ //ggrd_deal_with_composition_input(E, 0);
- } /* end grid cinit */
- /*
+ } /* end grid cinit */
+ /*
- check T and C range
+ check T and C range
- */
- if(E->control.check_t_irange)
- for(i=1;i<=E->lmesh.nno;i++){
- if(E->T[i]>1)E->T[i]=1;
- if(E->T[i]<0)E->T[i]=0;
- }
- if(E->control.check_c_irange)
- for(i=1;i<=E->lmesh.nno;i++){
- if(E->C[i]>1)E->C[i]=1;
- if(E->C[i]<0)E->C[i]=0;
- }
+ */
+ if(E->control.check_t_irange)
+ for(i=1;i<=E->lmesh.nno;i++){
+ if(E->T[i]>1)E->T[i]=1;
+ if(E->T[i]<0)E->T[i]=0;
+ }
+ if(E->control.check_c_irange)
+ for(i=1;i<=E->lmesh.nno;i++){
+ if(E->C[i]>1)E->C[i]=1;
+ if(E->C[i]<0)E->C[i]=0;
+ }
- if(E->control.composition)
- convection_initial_markers(E,1);
- } // end for restart==0
+ if(E->control.composition)
+ convection_initial_markers(E,1);
+ } // end for restart==0
else if(E->control.restart)
{
@@ -474,20 +476,20 @@
strncpy(ingfile,E->control.ggrd_flavor_gfile,1000);
strncpy(indfile,E->control.ggrd_flavor_dfile,1000);
if(E->tracers_add_flavors != 1)
- myerror("ggrd_deal_with_composition_input: only set up for one flavor",E);
+ myerror("ggrd_deal_with_composition_input: only set up for one flavor",E);
}
if(E->parallel.me > 0){
mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
MPI_COMM_WORLD, &mpi_stat);
}
- if(E->control.ggrd_slab_slice){
- if(E->control.ggrd_slab_slice == 1){
+ if(E->control.ggrd_c_slab_slice){
+ if(E->control.ggrd_c_slab_slice == 1){
if(ggrd_grdtrack_init_general(FALSE,ingfile, char_dummy,"",E->control.ggrd_ss_grd,(E->parallel.me==0),
FALSE,use_nearneighbor))
myerror("grdtrack x-z init error",E);
}else{ /* several slab slices */
- for(slice=0;slice < E->control.ggrd_slab_slice;slice++){
+ for(slice=0;slice < E->control.ggrd_c_slab_slice;slice++){
sprintf(pfile,"%s.%i.grd",ingfile,slice+1);
if(ggrd_grdtrack_init_general(FALSE,pfile, char_dummy,"",(E->control.ggrd_ss_grd+slice),(E->parallel.me==0),
FALSE,use_nearneighbor))
@@ -512,11 +514,11 @@
for(j=1;j<=nox;j++){
for(k=1;k<=noz;k++) {
node=k+(j-1)*noz+(i-1)*nox*noz; /* offset */
- if(E->control.ggrd_slab_slice){
+ if(E->control.ggrd_c_slab_slice){
/* slab */
- for(hit = slice=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
+ for(hit = slice=0;(!hit) && (slice < E->control.ggrd_c_slab_slice);slice++){
if(E->control.Rsphere) {
- if(in_slab_slice(E->SX[1][node],slice,E)){
+ if(in_slab_slice(E->SX[1][node],slice,E, E->control.ggrd_c_slab_theta_bound,E->control.ggrd_c_slab_slice)){
/* spherical interpolation for a slice */
ggrd_grdtrack_interpolate_xy((double)E->SX[2][node] * ONEEIGHTYOVERPI,
(double)E->SX[1][node],
@@ -525,7 +527,7 @@
hit = 1;
}
}else{ /* cartesian interpolation */
- if(in_slab_slice(E->X[2][node],slice,E)){
+ if(in_slab_slice(E->X[2][node],slice,E, E->control.ggrd_c_slab_theta_bound,E->control.ggrd_c_slab_slice)){
ggrd_grdtrack_interpolate_xy((double)E->X[1][node],(double)E->X[3][node],
(E->control.ggrd_ss_grd+slice),&tadd,
@@ -601,13 +603,13 @@
mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1), mpi_tag,
MPI_COMM_WORLD, &mpi_stat);
}
- if(E->control.ggrd_slab_slice){
- if(E->control.ggrd_slab_slice == 1){
+ if(E->control.ggrd_c_slab_slice){
+ if(E->control.ggrd_c_slab_slice == 1){
if(ggrd_grdtrack_init_general(FALSE,E->control.ggrd_flavor_gfile, char_dummy,"",
grd,(E->parallel.me==0),FALSE,use_nearneighbor))
myerror("grdtrack x-z init error",E);
}else{ /* several slab slices */
- for(slice=0;slice < E->control.ggrd_slab_slice;slice++){
+ for(slice=0;slice < E->control.ggrd_c_slab_slice;slice++){
sprintf(pfile,"%s.%i.grd",E->control.ggrd_flavor_gfile,slice+1);
if(ggrd_grdtrack_init_general(FALSE,pfile, char_dummy,"",(grd+slice),(E->parallel.me==0),
FALSE,use_nearneighbor))
@@ -630,11 +632,11 @@
fprintf(stderr,"ggrd_deal_with_composition_input: last processor done comp with grd init\n");
}
for(i=0;i < E->advection.markers;i++) {
- if(E->control.ggrd_slab_slice){
+ if(E->control.ggrd_c_slab_slice){
/* slab */
- for(hit = slice=0;(!hit) && (slice < E->control.ggrd_slab_slice);slice++){
+ for(hit = slice=0;(!hit) && (slice < E->control.ggrd_c_slab_slice);slice++){
if(E->control.Rsphere) {
- if(in_slab_slice(E->XMC[1][i],slice,E)){
+ if(in_slab_slice(E->XMC[1][i],slice,E, E->control.ggrd_c_slab_theta_bound,E->control.ggrd_c_slab_slice)){
/* spherical interpolation for a slice */
ggrd_grdtrack_interpolate_xy((double)E->XMC[2][i]* ONEEIGHTYOVERPI,(double)E->XMC[1][i],
(grd+slice),&tadd,
@@ -642,7 +644,7 @@
hit = 1;
}
}else{ /* cartesian interpolation */
- if(in_slab_slice(E->XMC[2][i],slice,E)){
+ if(in_slab_slice(E->XMC[2][i],slice,E, E->control.ggrd_c_slab_theta_bound,E->control.ggrd_c_slab_slice)){
ggrd_grdtrack_interpolate_xy((double)E->XMC[1][i],(double)E->XMC[3][i],
(grd+slice),&tadd,FALSE);
hit = 1 ;
@@ -681,26 +683,26 @@
-int in_slab_slice(float coord, int slice, struct All_variables *E)
+int in_slab_slice(float coord, int slice, struct All_variables *E,float *theta_bound, int slab_slice)
{
- if((slice < 0)||(slice > E->control.ggrd_slab_slice-1))
+ if((slice < 0)||(slice > slab_slice-1))
myerror("slab slice out of bounds",E);
- if(E->control.ggrd_slab_slice < 1)
+ if(slab_slice < 1)
myerror("total slab slice out of bounds",E);
- if(E->control.ggrd_slab_slice == 1)
- if(coord <= E->control.ggrd_slab_theta_bound[0])
+ if(slab_slice == 1)
+ if(coord <= theta_bound[0])
return 1;
else
return 0;
else{
if(slice == 0)
- if(coord <= E->control.ggrd_slab_theta_bound[0])
+ if(coord <= theta_bound[0])
return 1;
else
return 0;
else
- if((coord <= E->control.ggrd_slab_theta_bound[slice]) && (coord > E->control.ggrd_slab_theta_bound[slice-1]))
+ if((coord <= theta_bound[slice]) && (coord > theta_bound[slice-1]))
return 1;
else
return 0;
@@ -747,11 +749,11 @@
elxlylz = elxlz * ely;
/*
- if we have not initialized the time history structure, do it now
+ if we have not initialized the time history structure, do it now
*/
if(!E->control.ggrd.time_hist.init){
/*
- init times, if available
+ init times, if available
*/
ggrd_init_thist_from_file(&E->control.ggrd.time_hist,
E->control.ggrd.time_hist.file,TRUE,(E->parallel.me == 0));
@@ -770,7 +772,7 @@
mpi_rc = MPI_Recv(&mpi_inmsg, 1, MPI_INT, (E->parallel.me-1),
0, MPI_COMM_WORLD, &mpi_stat);
/*
- read in the material file(s)
+ read in the material file(s)
*/
E->control.ggrd.mat = (struct ggrd_gt *)calloc(E->control.ggrd.time_hist.nvtimes,sizeof(struct ggrd_gt));
/*
@@ -827,14 +829,14 @@
i1 = 0;
}
/*
- loop through all elements and assign
+ loop through all elements and assign
*/
for (j=1;j <= elz;j++) { /* this assumes a regular grid sorted as in (1)!!! */
if(((E->control.ggrd.mat_control > 0) && (E->mat[j] <= E->control.ggrd.mat_control )) ||
((E->control.ggrd.mat_control < 0) && (E->mat[j] == -E->control.ggrd.mat_control ))){
/*
lithosphere or asthenosphere
- */
+ */
for (k=1;k <= ely;k++){
for (i=1;i <= elx;i++) {
/* eq.(1) */
@@ -977,7 +979,7 @@
d[0] > d[1] > d[2]
- */
+*/
void ggrd_solve_eigen3x3(double a[3][3],double d[3],
double e[3][3],struct All_variables *E)
{
@@ -1104,7 +1106,7 @@
ntheta_grd,(E->parallel.me == 0),FALSE,FALSE))
myerror("ggrd init error",E);
if(E->control.CART3D)
- /* n_phi */
+ /* n_phi */
sprintf(tfilename,"%s/nx.grd",E->viscosity.anisotropic_init_dir);
else
sprintf(tfilename,"%s/np.grd",E->viscosity.anisotropic_init_dir);
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2011-06-11 00:25:45 UTC (rev 18590)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2011-06-12 02:11:10 UTC (rev 18591)
@@ -678,6 +678,7 @@
/* first the problem type (defines subsequent behaviour) */
int m;
+ int i;
m = E->parallel.me;
@@ -857,11 +858,39 @@
input_double("ggrd_cinit_offset",&(E->control.ggrd.comp.offset),"0.0", m);
- /* slab slice handling */
- input_int("slab_slice",&(E->control.ggrd_slab_slice),"0", m);
- if(E->control.ggrd_slab_slice > 3)
- myerror("too many slab slices, only four allowed",E);
- input_float_vector("slab_theta_bound",E->control.ggrd_slab_slice,(E->control.ggrd_slab_theta_bound), m);
+ /*
+
+ slab slice handling
+
+ */
+ input_int("slab_slice",&(E->control.ggrd_t_slab_slice),"0", m);
+ /* old mode, both composition and temperature are controlled
+ by the same number of slices and theta bounds */
+ if(E->control.ggrd_t_slab_slice){
+ if(m == 0)fprintf(stderr,"WARNING: compatibility mode, both comp and temp use the same slice number and geometry\n");
+ if(E->control.ggrd_t_slab_slice > GGRD_MAX_NR_SLICE)
+ myerror("too many slab slices, increase GGRD_MAX_NR_SLICE",E);
+ E->control.ggrd_c_slab_slice = E->control.ggrd_t_slab_slice;
+ input_float_vector("slab_theta_bound",E->control.ggrd_t_slab_slice,(E->control.ggrd_t_slab_theta_bound), m);
+ for(i=0;i<GGRD_MAX_NR_SLICE;i++){
+ E->control.ggrd_c_slab_theta_bound[i] = E->control.ggrd_t_slab_theta_bound[i];
+ }
+ }else{
+ /* allow for different slices */
+ input_int("slab_t_slice",&(E->control.ggrd_t_slab_slice),"0", m);
+ if(E->control.ggrd_t_slab_slice > GGRD_MAX_NR_SLICE)
+ myerror("too many temp slab slices, increase GGRD_MAX_NR_SLICE",E);
+ input_float_vector("slab_t_theta_bound",E->control.ggrd_t_slab_slice,(E->control.ggrd_t_slab_theta_bound), m);
+
+ input_int("slab_c_slice",&(E->control.ggrd_c_slab_slice),"0", m);
+ if(E->control.ggrd_c_slab_slice > GGRD_MAX_NR_SLICE)
+ myerror("too many comp slab slices, increase GGRD_MAX_NR_SLICE",E);
+ input_float_vector("slab_c_theta_bound",E->control.ggrd_c_slab_slice,(E->control.ggrd_c_slab_theta_bound), m);
+
+
+
+
+ }
/*
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2011-06-11 00:25:45 UTC (rev 18590)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2011-06-12 02:11:10 UTC (rev 18591)
@@ -94,9 +94,11 @@
#define MAX_NEIGHBORS 27
-//#define CU_MPI_MSG_LIM 100
-#define CU_MPI_MSG_LIM 1000
+#define GGRD_MAX_NR_SLICE 5
+#define CU_MPI_MSG_LIM 100 /* this increase wasn't necessary */
+//#define CU_MPI_MSG_LIM 1000
+
/* Macros */
#define max(A,B) (((A) > (B)) ? (A) : (B))
@@ -763,8 +765,10 @@
float Ra_410, clapeyron410, transT410, width410;
#ifdef USE_GGRD
struct ggrd_master ggrd;
- int ggrd_slab_slice;
- float ggrd_slab_theta_bound[4];
+
+ int ggrd_t_slab_slice,ggrd_c_slab_slice;
+ float ggrd_t_slab_theta_bound[GGRD_MAX_NR_SLICE],ggrd_c_slab_theta_bound[GGRD_MAX_NR_SLICE];
+
struct ggrd_gt ggrd_ss_grd[4];
int ggrd_mat_limit_prefactor,ggrd_mat_is_3d;
char ggrd_mat_depth_file[1000],ggrd_flavor_gfile[1000],ggrd_flavor_dfile[1000];
@@ -1063,11 +1067,14 @@
could also use
cproto -DCITCOM_ALLOW_ANISOTROPIC_VISC -DUSE_GGRD -q -p -I. -f2 *.c > prototypes.h
- */
+ cproto -DCITCOM_ALLOW_ANISOTROPIC_VISC -DUSE_GGRD -q -p -I. -I$HOME/progs/src/hc-svn/ -I$GMTHOME/include/ -f2 *.c > prototypes.h
+
+
+ */
void twiddle_thumbs(struct All_variables *, int ); /* somehow, cproto
does not pick up
these functions?! */
-void convection_initial_temperature_ggrd(struct All_variables *);
+void convection_initial_temperature_and_comp_ggrd(struct All_variables *);
void set_2dc_defaults(struct All_variables *);
void remove_horiz_ave(struct All_variables *, float *, float *, int );
void output_velo_related_gzdir(struct All_variables *, int);
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2011-06-11 00:25:45 UTC (rev 18590)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2011-06-12 02:11:10 UTC (rev 18591)
@@ -39,9 +39,10 @@
void f_times_ft(double [3][3], double [3][3]);
void drex_eigen(double [3][3], double [3][3], int *);
void malmul_scaled_id(double [3][3], double [3][3], double, double);
+void myerror_s(char *, struct All_variables *);
/* Boundary_conditions.c */
+void velocity_boundary_conditions(struct All_variables *);
void freeze_surface(struct All_variables *);
-void velocity_boundary_conditions(struct All_variables *);
void temperature_boundary_conditions(struct All_variables *);
void velocity_refl_vert_bc(struct All_variables *);
void temperature_refl_vert_bc(struct All_variables *);
@@ -63,6 +64,7 @@
void prepare_transfer_arrays(struct All_variables *);
int locate_processor(struct All_variables *, double, double, double);
void get_C_from_markers(struct All_variables *, float *);
+void get_CF_from_markers(struct All_variables *, int **);
void element_markers(struct All_variables *, int);
void velocity_markers(struct All_variables *, float *[4], int);
int get_element(struct All_variables *, double, double, double, double [4]);
@@ -97,6 +99,7 @@
void process_restart_tc(struct All_variables *);
void convection_initial_markers1(struct All_variables *);
void convection_initial_markers(struct All_variables *, int);
+void assign_flavor_to_tracer_based_on_node(struct All_variables *, int, int);
void setup_plume_problem(struct All_variables *);
void PG_process(struct All_variables *, int);
/* Drive_solvers.c */
@@ -107,7 +110,6 @@
void get_elt_k(struct All_variables *, int, double [24 * 24], int, int);
void get_ba(struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [9][4][9][7]);
void get_ba_p(struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [9][4][9][7]);
-void get_vgm_p(double [4][9], struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [3][3], double [3]);
void assemble_del2_u(struct All_variables *, double *, double *, int, int);
void e_assemble_del2_u(struct All_variables *, double *, double *, int, int);
void n_assemble_del2_u(struct All_variables *, double *, double *, int, int);
@@ -166,12 +168,14 @@
void set_3ds_defaults(struct All_variables *);
void set_3dc_defaults(struct All_variables *);
/* Ggrd_handling.c */
-int in_slab_slice(float, int, struct All_variables *);
+void convection_initial_temperature_and_comp_ggrd(struct All_variables *);
+void assign_flavor_to_tracer_from_grd(struct All_variables *);
+void ggrd_deal_with_composition_input(struct All_variables *, int);
+void assign_flavor_to_tracer_from_grd(struct All_variables *);
+int in_slab_slice(float, int, struct All_variables *,float *,int);
void ggrd_read_mat_from_file(struct All_variables *);
void ggrd_solve_eigen3x3(double [3][3], double [3], double [3][3], struct All_variables *);
void ggrd_read_anivisc_from_file(struct All_variables *);
-void ggrd_deal_with_composition_input(struct All_variables *, int );
-
/* Global_operations.c */
void return_horiz_sum(struct All_variables *, float *, float *, int);
void return_horiz_ave(struct All_variables *, float *, float *);
@@ -203,14 +207,20 @@
void common_initial_fields(struct All_variables *);
void initial_pressure(struct All_variables *);
void initial_velocity(struct All_variables *);
+/* io.c */
+_Bool io_directory_create(const char *);
+_Bool io_directory_exists(const char *);
+void io_error(void);
+_Bool io_results(struct All_variables *);
+_Bool io_setup_path(struct All_variables *);
/* Nodal_mesh.c */
void node_locations(struct All_variables *);
void pre_interpolation(struct All_variables *);
void dlogical_mesh_to_real(struct All_variables *, double *, int);
void flogical_mesh_to_real(struct All_variables *, float *, int);
void p_to_nodes(struct All_variables *, double *, float *, int);
+void e2_to_nodes(struct All_variables *, float *, float *, int);
void p_to_centres(struct All_variables *, float *, double *, int);
-void e2_to_nodes(struct All_variables *, float *, float *, int );
void v_to_intpts(struct All_variables *, float *, float *, int);
void v_to_nodes(struct All_variables *, float *, float *, int);
void visc_to_intpts(struct All_variables *, float *, float *, int);
@@ -219,6 +229,8 @@
void visc_from_gint_to_ele(struct All_variables *, float *, float *, int);
void visc_from_gint_to_nodes(struct All_variables *, float *, float *, int);
void visc_from_nodes_to_gint(struct All_variables *, float *, float *, int);
+/* output_ascii.c */
+_Bool output_ascii(struct All_variables *);
/* Output.c */
void output_velo_related_binary(struct All_variables *, int);
void output_temp(struct All_variables *, int);
@@ -227,6 +239,8 @@
/* Output_gzdir.c */
void output_velo_related_gzdir(struct All_variables *, int);
void process_restart_tc_gzdir(struct All_variables *);
+/* output_vtk.c */
+_Bool output_vtk(struct All_variables *);
/* Pan_problem_misc_functions.c */
int get_process_identifier(void);
void unique_copy_file(struct All_variables *, char *, char *);
@@ -261,6 +275,7 @@
void matmul_3x3(double [3][3], double [3][3], double [3][3]);
void assign_to_3x3(double [3][3], double);
void remove_trace_3x3(double [3][3]);
+double distance_to_node(double, double, double, struct All_variables *, int);
/* Parallel_related.c */
void parallel_process_termination(void);
void parallel_domain_decomp1(struct All_variables *);
@@ -276,6 +291,7 @@
void exchange_markers(struct All_variables *);
void exchange_id_d20(struct All_variables *, double *, int);
void exchange_node_f20(struct All_variables *, float *, int);
+void exchange_node_int(struct All_variables *, int *, int);
double CPU_time0(void);
void parallel_process_sync(void);
/* Parsing.c */
@@ -294,14 +310,11 @@
int input_double_vector(char *, int, double *, int);
int interpret_control_string(char *, int *, double *, double *, double *);
/* Phase_change.c */
-void phase_change(struct All_variables *, float *, float *, float *, float *);
/* Process_buoyancy.c */
-void process_temp_field(struct All_variables *, int);
void heat_flux(struct All_variables *);
void heat_flux1(struct All_variables *);
void plume_buoyancy_flux(struct All_variables *);
/* Process_velocity.c */
-void process_new_velocity(struct All_variables *, int);
void get_surface_velo(struct All_variables *, float *);
void get_ele_visc(struct All_variables *, float *);
void get_surf_stress(struct All_variables *, float *, float *, float *, float *, float *, float *);
@@ -309,11 +322,9 @@
/* Profiling.c */
float CPU_time(void);
/* Shape_functions.c */
-void construct_shape_functions(struct All_variables *);
double lpoly(int, double);
double lpolydash(int, double);
/* Size_does_matter.c */
-void twiddle_thumbs(struct All_variables *, int);
void get_global_shape_fn(struct All_variables *, int, int, double [4][9], int, int);
void form_rtf_bc(int, double [4], double [4][9], double [4][4]);
void get_rtf(struct All_variables *, int, int, double [4][9], int);
@@ -322,7 +333,6 @@
void get_global_1d_shape_fn1(struct All_variables *, int, struct Shape_function1 *, struct Shape_function1_dA *, int);
void mass_matrix(struct All_variables *);
/* Solver_conj_grad.c */
-void set_cg_defaults(struct All_variables *);
void cg_allocate_vars(struct All_variables *);
void assemble_forces_iterative(struct All_variables *);
/* Solver_multigrid.c */
@@ -339,7 +349,6 @@
void inject_scalar(struct All_variables *, int, float *, float *);
void inject_scalar_e(struct All_variables *, int, float *, float *);
/* Sphere_harmonics.c */
-void set_sphere_harmonics(struct All_variables *);
void sphere_harmonics_layer(struct All_variables *, float **, float *, float *, int, char *);
void sphere_interpolate(struct All_variables *, float **, float *);
void sphere_expansion(struct All_variables *, float *, float *, float *);
@@ -364,10 +373,10 @@
void visc_from_mat(struct All_variables *, float *, float *);
void visc_from_T(struct All_variables *, float *, float *, int);
void visc_from_S(struct All_variables *, float *, float *, int);
-void visc_from_S2(struct All_variables *, float *, float *, int);
void strain_rate_2_inv(struct All_variables *, float *, int);
double second_invariant_from_3x3(double [3][3]);
void calc_vgm_matrix(struct All_variables *, double *, double *);
+void get_vgm_p(double [4][9], struct Shape_function *, struct Shape_function_dx *, struct CC *, struct CCX *, double [4][9], int, int, int, int, double [3][3], double [3]);
void calc_strain_from_vgm(double [3][3], double [3][3]);
void calc_strain_from_vgm9(double *, double [3][3]);
void calc_rot_from_vgm(double [3][3], double [3][3]);
@@ -375,8 +384,3 @@
int layers(struct All_variables *, float);
int weak_zones(struct All_variables *, int, float);
float boundary_thickness(struct All_variables *, float *);
-double distance_to_node(double , double , double ,struct All_variables *, int );
-void assign_flavor_to_tracer_based_on_node(struct All_variables *, int , int );
-void exchange_node_int(struct All_variables *, int *, int );
-void get_CF_from_markers(struct All_variables *, int **);
-void assign_flavor_to_tracer_from_grd(struct All_variables *);
More information about the CIG-COMMITS
mailing list