[cig-commits] r18754 - in mc/3D/CitcomS/branches/eheien: . bin lib module visual
emheien at geodynamics.org
emheien at geodynamics.org
Wed Jul 13 17:06:31 PDT 2011
Author: emheien
Date: 2011-07-13 17:06:30 -0700 (Wed, 13 Jul 2011)
New Revision: 18754
Modified:
mc/3D/CitcomS/branches/eheien/bin/Citcom.c
mc/3D/CitcomS/branches/eheien/bin/CitcomSFull.c
mc/3D/CitcomS/branches/eheien/bin/CitcomSRegional.c
mc/3D/CitcomS/branches/eheien/configure.ac
mc/3D/CitcomS/branches/eheien/lib/Advection_diffusion.c
mc/3D/CitcomS/branches/eheien/lib/BC_util.c
mc/3D/CitcomS/branches/eheien/lib/Checkpoints.c
mc/3D/CitcomS/branches/eheien/lib/Composition_related.c
mc/3D/CitcomS/branches/eheien/lib/Construct_arrays.c
mc/3D/CitcomS/branches/eheien/lib/Convection.c
mc/3D/CitcomS/branches/eheien/lib/Element_calculations.c
mc/3D/CitcomS/branches/eheien/lib/Full_boundary_conditions.c
mc/3D/CitcomS/branches/eheien/lib/Full_geometry_cartesian.c
mc/3D/CitcomS/branches/eheien/lib/Full_parallel_related.c
mc/3D/CitcomS/branches/eheien/lib/Full_read_input_from_files.c
mc/3D/CitcomS/branches/eheien/lib/Full_sphere_related.c
mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c
mc/3D/CitcomS/branches/eheien/lib/Full_version_dependent.c
mc/3D/CitcomS/branches/eheien/lib/General_matrix_functions.c
mc/3D/CitcomS/branches/eheien/lib/Ggrd_handling.c
mc/3D/CitcomS/branches/eheien/lib/Global_operations.c
mc/3D/CitcomS/branches/eheien/lib/Initial_temperature.c
mc/3D/CitcomS/branches/eheien/lib/Instructions.c
mc/3D/CitcomS/branches/eheien/lib/Interuption.c
mc/3D/CitcomS/branches/eheien/lib/Makefile.am
mc/3D/CitcomS/branches/eheien/lib/Mineral_physics_models.c
mc/3D/CitcomS/branches/eheien/lib/Nodal_mesh.c
mc/3D/CitcomS/branches/eheien/lib/Output.c
mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c
mc/3D/CitcomS/branches/eheien/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/eheien/lib/Parsing.c
mc/3D/CitcomS/branches/eheien/lib/Problem_related.c
mc/3D/CitcomS/branches/eheien/lib/Process_buoyancy.c
mc/3D/CitcomS/branches/eheien/lib/Regional_boundary_conditions.c
mc/3D/CitcomS/branches/eheien/lib/Regional_geometry_cartesian.c
mc/3D/CitcomS/branches/eheien/lib/Regional_parallel_related.c
mc/3D/CitcomS/branches/eheien/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/branches/eheien/lib/Regional_sphere_related.c
mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/eheien/lib/Regional_version_dependent.c
mc/3D/CitcomS/branches/eheien/lib/Shape_functions.c
mc/3D/CitcomS/branches/eheien/lib/Size_does_matter.c
mc/3D/CitcomS/branches/eheien/lib/Solver_conj_grad.c
mc/3D/CitcomS/branches/eheien/lib/Solver_multigrid.c
mc/3D/CitcomS/branches/eheien/lib/Sphere_harmonics.c
mc/3D/CitcomS/branches/eheien/lib/Sphere_util.c
mc/3D/CitcomS/branches/eheien/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/branches/eheien/lib/Topo_gravity.c
mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c
mc/3D/CitcomS/branches/eheien/lib/global_defs.h
mc/3D/CitcomS/branches/eheien/lib/prototypes.h
mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h
mc/3D/CitcomS/branches/eheien/module/misc.c
mc/3D/CitcomS/branches/eheien/visual/project_geoid.c
Log:
Initial work in converting tracer code to C++
Changed autoconf to use C++ compiler
Fixed numerous C++ syntax violations
Fixed several minor bugs
Modified: mc/3D/CitcomS/branches/eheien/bin/Citcom.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/bin/Citcom.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/bin/Citcom.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -42,9 +42,7 @@
void solver_init(struct All_variables *E);
-int main(argc,argv)
- int argc;
- char **argv;
+int main(int argc, char **argv)
{ /* Functions called by main*/
void general_stokes_solver();
Modified: mc/3D/CitcomS/branches/eheien/bin/CitcomSFull.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/bin/CitcomSFull.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/bin/CitcomSFull.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -28,8 +28,16 @@
struct All_variables;
+#ifdef __cplusplus
+extern "C" {
+#endif
+
void full_solver_init(struct All_variables *E);
+#ifdef __cplusplus
+}
+#endif
+
void solver_init(struct All_variables *E)
{
full_solver_init(E);
Modified: mc/3D/CitcomS/branches/eheien/bin/CitcomSRegional.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/bin/CitcomSRegional.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/bin/CitcomSRegional.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -28,8 +28,16 @@
struct All_variables;
+#ifdef __cplusplus
+extern "C" {
+#endif
+
void regional_solver_init(struct All_variables *E);
+#ifdef __cplusplus
+}
+#endif
+
void solver_init(struct All_variables *E)
{
regional_solver_init(E);
Modified: mc/3D/CitcomS/branches/eheien/configure.ac
===================================================================
--- mc/3D/CitcomS/branches/eheien/configure.ac 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/configure.ac 2011-07-14 00:06:30 UTC (rev 18754)
@@ -71,6 +71,9 @@
AC_PROG_CC([mpicc hcc mpcc mpcc_r mpxlc cmpicc gcc cc cl icc ecc pgcc xlc xlc_r])
AC_PROG_CXX([mpicxx mpic++ mpiCC hcp mpCC mpxlC mpxlC_r cmpic++ g++ c++ gpp aCC CC cxx cc++ cl FCC KCC RCC xlc++_r xlC_r xlC icpc ecpc pgCC])
+CC="$CXX"
+CFLAGS="$CXXFLAGS"
+
# Checking a system header here so that CPP is always defined
AC_CHECK_HEADERS([malloc.h])
Modified: mc/3D/CitcomS/branches/eheien/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Advection_diffusion.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Advection_diffusion.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -143,7 +143,6 @@
{
int i,d,n,nel,el,node,m;
- float global_fmin();
void velo_from_element();
float adv_timestep;
@@ -178,12 +177,12 @@
uc = fabs(uc1)/E->eco[m][el].size[1] + fabs(uc2)/E->eco[m][el].size[2] + fabs(uc3)/E->eco[m][el].size[3];
step = (0.5/uc);
- adv_timestep = min(adv_timestep,step);
+ adv_timestep = fmin(adv_timestep,step);
}
adv_timestep = E->advection.dt_reduced * adv_timestep;
- adv_timestep = 1.0e-32 + min(E->advection.fine_tune_dt*adv_timestep,
+ adv_timestep = 1.0e-32 + fmin(E->advection.fine_tune_dt*adv_timestep,
E->advection.diff_timestep);
E->advection.timestep = global_fmin(E,adv_timestep);
@@ -246,7 +245,7 @@
process_heating(E, psc_pass);
/* XXX: replace inputdiff with refstate.thermal_conductivity */
- pg_solver(E,E->T,E->Tdot,DTdot,&(E->convection.heat_sources),E->control.inputdiff,1,E->node);
+ pg_solver(E,E->T,E->Tdot,DTdot,(SOURCES*)&(E->convection.heat_sources),E->control.inputdiff,1,E->node);
corrector(E,E->T,E->Tdot,DTdot);
temperatures_conform_bcs(E);
}
@@ -320,14 +319,12 @@
float diff_timestep, ts;
int m, el, d;
- float global_fmin();
-
diff_timestep = 1.0e8;
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(el=1;el<=E->lmesh.nel;el++) {
for(d=1;d<=E->mesh.nsd;d++) {
ts = E->eco[m][el].size[d] * E->eco[m][el].size[d];
- diff_timestep = min(diff_timestep,ts);
+ diff_timestep = fmin(diff_timestep,ts);
}
}
Modified: mc/3D/CitcomS/branches/eheien/lib/BC_util.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/BC_util.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/BC_util.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -33,6 +33,8 @@
int layers(struct All_variables *,int,int);
+#include <math.h>
+
#ifdef USE_GGRD
#include "ggrd_handling.h"
#endif
@@ -86,10 +88,7 @@
}
-void strip_bcs_from_residual(E,Res,level)
- struct All_variables *E;
- double **Res;
- int level;
+void strip_bcs_from_residual(struct All_variables *E, double **Res, int level)
{
int m,i;
@@ -102,8 +101,7 @@
}
-void temperatures_conform_bcs(E)
- struct All_variables *E;
+void temperatures_conform_bcs(struct All_variables *E)
{
void temperatures_conform_bcs2(struct All_variables *);
void assimilate_lith_conform_bcs2(struct All_variables *);
@@ -121,8 +119,7 @@
}
-void temperatures_conform_bcs2(E)
- struct All_variables *E;
+void temperatures_conform_bcs2(struct All_variables *E)
{
int j,node;
unsigned int type;
@@ -166,9 +163,7 @@
}
-void velocities_conform_bcs(E,U)
- struct All_variables *E;
- double **U;
+void velocities_conform_bcs(struct All_variables *E, double **U)
{
int node,m;
Modified: mc/3D/CitcomS/branches/eheien/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Checkpoints.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Checkpoints.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -251,6 +251,7 @@
write_sentinel(fp);
+ /*
fwrite(&(E->trace.number_of_basic_quantities), sizeof(int), 1, fp);
fwrite(&(E->trace.number_of_extra_quantities), sizeof(int), 1, fp);
fwrite(&(E->trace.nflavors), sizeof(int), 1, fp);
@@ -258,9 +259,11 @@
for(m=1; m<=E->sphere.caps_per_proc; m++)
fwrite(&(E->trace.ntracers[m]), sizeof(int), 1, fp);
+ */
/* the 0-th element of basicq/extraq/ielement is not init'd
* and won't be used when read it. */
+ /*
for(m=1; m<=E->sphere.caps_per_proc; m++) {
for(i=0; i<6; i++) {
fwrite(E->trace.basicq[m][i], sizeof(double),
@@ -273,6 +276,7 @@
fwrite(E->trace.ielement[m], sizeof(int),
E->trace.ntracers[m]+1, fp);
}
+ */
return;
}
@@ -287,6 +291,7 @@
read_sentinel(fp, E->parallel.me);
+ /*
fread(&itmp, sizeof(int), 1, fp);
if (itmp != E->trace.number_of_basic_quantities) {
fprintf(stderr, "Error in reading checkpoint file: tracer basicq, me=%d\n",
@@ -317,14 +322,19 @@
fread(&itmp, sizeof(int), 1, fp);
E->trace.ilast_tracer_count = itmp;
+ */
+
/* # of tracers, allocate memory */
+ /*
for(m=1; m<=E->sphere.caps_per_proc; m++) {
fread(&itmp, sizeof(int), 1, fp);
allocate_tracer_arrays(E, m, itmp);
E->trace.ntracers[m] = itmp;
}
+ */
/* read tracer data */
+ /*
for(m=1; m<=E->sphere.caps_per_proc; m++) {
for(i=0; i<6; i++) {
fread(E->trace.basicq[m][i], sizeof(double),
@@ -337,6 +347,7 @@
fread(E->trace.ielement[m], sizeof(int),
E->trace.ntracers[m]+1, fp);
}
+ */
/* init E->trace.ntracer_flavor */
count_tracers_of_flavors(E);
Modified: mc/3D/CitcomS/branches/eheien/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Composition_related.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Composition_related.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -303,7 +303,7 @@
for (e=1; e<=E->lmesh.nel; e++) {
numtracers = 0;
for (flavor=0; flavor<E->trace.nflavors; flavor++)
- numtracers += E->trace.ntracer_flavor[j][flavor][e];
+ numtracers += E->trace.num_tracer_flavors[j][flavor][e];
/* Check for empty entries and compute ratio. */
/* If no tracers are in an element, skip this element, */
@@ -317,7 +317,7 @@
for(i=0;i<E->composition.ncomp;i++) {
flavor = i + 1;
E->composition.comp_el[j][i][e] =
- E->trace.ntracer_flavor[j][flavor][e] / (double)numtracers;
+ E->trace.num_tracer_flavors[j][flavor][e] / (double)numtracers;
}
}
@@ -370,7 +370,7 @@
for (e=1; e<=E->lmesh.nel; e++) {
numtracers = 0;
for (flavor=0; flavor<E->trace.nflavors; flavor++)
- numtracers += E->trace.ntracer_flavor[j][flavor][e];
+ numtracers += E->trace.num_tracer_flavors[j][flavor][e];
/* Check for empty entries */
/* If no tracers are in an element, comp = 0.0 (i.e. is ambient) */
@@ -385,13 +385,13 @@
for(i=0;i<E->composition.ncomp;i++) {
flavor = i;
comp =
- E->trace.ntracer_flavor[j][flavor][e] / E->eco[j][e].area
- * domain_volume / E->trace.number_of_tracers;
+ E->trace.num_tracer_flavors[j][flavor][e] / E->eco[j][e].area
+ * domain_volume / E->trace.tracers.size();
/* truncate composition at 1.0 */
/* This violates mass conservation but prevents unphysical C */
/* XXX: make truncation a switch for the user to specify */
- E->composition.comp_el[j][i][e] = min(comp,one);
+ E->composition.comp_el[j][i][e] = fmin(comp,one);
}
}
@@ -503,7 +503,7 @@
for (e=1; e<=E->lmesh.nel; e++) {
numtracers = 0;
for (flavor=0; flavor<E->trace.nflavors; flavor++)
- numtracers += E->trace.ntracer_flavor[j][flavor][e];
+ numtracers += E->trace.num_tracer_flavors[j][flavor][e];
if (numtracers == 0)
is_empty[e] = 1;
Modified: mc/3D/CitcomS/branches/eheien/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Construct_arrays.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Construct_arrays.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -42,9 +42,7 @@
it should be done through a pre-calculated lookup table.
======================================================== */
-void construct_ien(E)
- struct All_variables *E;
-
+void construct_ien(struct All_variables *E)
{
int lev,p,q,r,rr,j;
int element,start,nel,nno;
@@ -136,8 +134,7 @@
Function to make the ID array for above case
============================================ */
-void construct_id(E)
- struct All_variables *E;
+void construct_id(struct All_variables *E)
{
int i,j,k;
int eqn_count,node,nno;
@@ -199,9 +196,7 @@
-void get_bcs_id_for_residual(E,level,m)
- struct All_variables *E;
- int level,m;
+void get_bcs_id_for_residual(struct All_variables *E, int level, int m)
{
int i,j;
@@ -233,8 +228,7 @@
Function to construct the LM array from the ID and IEN arrays
========================================================== */
-void construct_lm(E)
- struct All_variables *E;
+void construct_lm(struct All_variables *E)
{
int i,j,a,e;
int lev,eqn_no;
@@ -251,8 +245,7 @@
Function to build the local node matrix indexing maps
===================================================== */
-void construct_node_maps(E)
- struct All_variables *E;
+void construct_node_maps(struct All_variables *E)
{
double time1,CPU_time0();
@@ -328,8 +321,7 @@
}
-void construct_node_ks(E)
- struct All_variables *E;
+void construct_node_ks(struct All_variables *E)
{
int m,level,i,j,k,e;
int node,node1,eqn1,eqn2,eqn3,loc0,loc1,loc2,loc3,found,element,index,pp,qq;
@@ -476,8 +468,7 @@
return;
}
-void rebuild_BI_on_boundary(E)
- struct All_variables *E;
+void rebuild_BI_on_boundary(struct All_variables *E)
{
int m,level,i,j;
int eqn1,eqn2,eqn3;
@@ -545,8 +536,7 @@
masks and other indicators.
============================================ */
-void construct_masks(E) /* Add lid/edge masks/nodal weightings */
- struct All_variables *E;
+void construct_masks(struct All_variables *E) /* Add lid/edge masks/nodal weightings */
{
int i,j,k,l,node,el,elt;
int lev,elx,elz,ely,nno,nox,noz,noy;
@@ -595,8 +585,7 @@
build the sub-element reference matrices
========================================== */
-void construct_sub_element(E)
- struct All_variables *E;
+void construct_sub_element(struct All_variables *E)
{ int i,j,k,l,m;
int lev,nox,noy,noz,nnn,elx,elz,ely,elzu,elxu,elt,eltu;
@@ -641,8 +630,7 @@
}
-void construct_elt_ks(E)
- struct All_variables *E;
+void construct_elt_ks(struct All_variables *E)
{
int e,el,lev,j,k,ii,m;
void get_elt_k();
@@ -688,8 +676,7 @@
-void construct_elt_gs(E)
- struct All_variables *E;
+void construct_elt_gs(struct All_variables *E)
{ int m,el,lev,a;
void get_elt_g();
@@ -737,8 +724,7 @@
routine for constructing stiffness and node_maps
============================================================== */
-void construct_stiffness_B_matrix(E)
- struct All_variables *E;
+void construct_stiffness_B_matrix(struct All_variables *E)
{
void build_diagonal_of_K();
void build_diagonal_of_Ahat();
@@ -818,8 +804,7 @@
============================================================== */
-void construct_mat_group(E)
- struct All_variables *E;
+void construct_mat_group(struct All_variables *E)
{
int m,i,j,k,kk,el,lev,a,nodea,els,llayer;
void read_visc_layer_file(struct All_variables *E);
Modified: mc/3D/CitcomS/branches/eheien/lib/Convection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Convection.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Convection.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -31,8 +31,7 @@
#include "global_defs.h"
-void set_convection_defaults(E)
- struct All_variables *E;
+void set_convection_defaults(struct All_variables *E)
{
/*
void PG_timestep();
@@ -71,9 +70,7 @@
return;
}
-void read_convection_settings(E)
- struct All_variables *E;
-
+void read_convection_settings(struct All_variables *E)
{
void advection_diffusion_parameters();
@@ -88,16 +85,13 @@
Any setup which relates only to the convection stuff goes in here
================================================================= */
-void convection_derived_values(E)
- struct All_variables *E;
-
+void convection_derived_values(struct All_variables *E)
{
return;
}
-void convection_allocate_memory(E)
- struct All_variables *E;
+void convection_allocate_memory(struct All_variables *E)
{ void advection_diffusion_allocate_memory();
@@ -108,8 +102,7 @@
/* ============================================ */
-void convection_initial_fields(E)
- struct All_variables *E;
+void convection_initial_fields(struct All_variables *E)
{
void convection_initial_temperature();
@@ -120,8 +113,7 @@
/* =========================================== */
-void convection_boundary_conditions(E)
- struct All_variables *E;
+void convection_boundary_conditions(struct All_variables *E)
{
(E->solver.velocity_boundary_conditions)(E); /* universal */
Modified: mc/3D/CitcomS/branches/eheien/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Element_calculations.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Element_calculations.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -76,9 +76,7 @@
Function to get the global H vector (mixed method driving terms)
================================================================ */
-void assemble_forces(E,penalty)
- struct All_variables *E;
- int penalty;
+void assemble_forces(struct All_variables *E, int penalty)
{
double elt_f[24];
int m,a,e,i;
@@ -278,11 +276,7 @@
Function to supply the element k matrix for a given element e.
============================================================== */
-void get_elt_k(E,el,elt_k,lev,m,iconv)
- struct All_variables *E;
- int el,m;
- double elt_k[24*24];
- int lev, iconv;
+void get_elt_k(struct All_variables *E, int el, double elt_k[24*24], int lev, int m, int iconv)
{
double bdbmu[4][4];
int pn,qn,ad,bd,off;
@@ -424,11 +418,7 @@
element or node by node.
============================================= */
-void assemble_del2_u(E,u,Au,level,strip_bcs)
- struct All_variables *E;
- double **u,**Au;
- int level;
- int strip_bcs;
+void assemble_del2_u(struct All_variables *E, double **u, double **Au, int level, int strip_bcs)
{
void e_assemble_del2_u();
void n_assemble_del2_u();
@@ -445,11 +435,7 @@
Assemble del_squared_u vector el by el
====================================== */
-void e_assemble_del2_u(E,u,Au,level,strip_bcs)
- struct All_variables *E;
- double **u,**Au;
- int level;
- int strip_bcs;
+void e_assemble_del2_u(struct All_variables *E, double **u, double **Au, int level, int strip_bcs)
{
int e,i,a,b,a1,a2,a3,ii,m,nodeb;
@@ -518,11 +504,7 @@
Assemble Au using stored, nodal coefficients.
====================================================== */
-void n_assemble_del2_u(E,u,Au,level,strip_bcs)
- struct All_variables *E;
- double **u,**Au;
- int level;
- int strip_bcs;
+void n_assemble_del2_u(struct All_variables *E, double **u, double **Au, int level, int strip_bcs)
{
int m, e,i;
int eqn1,eqn2,eqn3;
@@ -582,11 +564,7 @@
}
-void build_diagonal_of_K(E,el,elt_k,level,m)
- struct All_variables *E;
- int level,el,m;
- double elt_k[24*24];
-
+void build_diagonal_of_K(struct All_variables *E, int el, double elt_k[24*24], int level, int m)
{
int a,a1,a2,p,node;
@@ -615,8 +593,7 @@
return;
}
-void build_diagonal_of_Ahat(E)
- struct All_variables *E;
+void build_diagonal_of_Ahat(struct All_variables *E)
{
double assemble_dAhatp_entry();
@@ -739,11 +716,7 @@
Assemble a grad_P vector element by element
========================================== */
-void assemble_grad_p(E,P,gradP,lev)
- struct All_variables *E;
- double **P,**gradP;
- int lev;
-
+void assemble_grad_p(struct All_variables *E, double **P, double **gradP, int lev)
{
int m,e,i,j1,j2,j3,p,a,b,nel,neq;
void strip_bcs_from_residual();
@@ -787,10 +760,7 @@
}
-double assemble_dAhatp_entry(E,e,level,m)
- struct All_variables *E;
- int e,level,m;
-
+double assemble_dAhatp_entry(struct All_variables *E, int e, int level, int m)
{
int i,j,p,a,b,node,npno;
void strip_bcs_from_residual();
@@ -915,12 +885,7 @@
used for the divergence
============================================================== */
-void get_elt_g(E,el,elt_del,lev,m)
- struct All_variables *E;
- int el,m;
- higher_precision elt_del[24][1];
- int lev;
-
+void get_elt_g(struct All_variables *E, int el, higher_precision elt_del[24][1], int lev, int m)
{
void get_rtf_at_ppts();
int p,a,i,j,k;
@@ -1017,12 +982,7 @@
Function to create the element force vector (allowing for velocity b.c.'s)
================================================================= */
-void get_elt_f(E,el,elt_f,bcs,m)
- struct All_variables *E;
- int el,m;
- double elt_f[24];
- int bcs;
-
+void get_elt_f(struct All_variables *E, int el, double elt_f[24], int bcs, int m)
{
int i,p,a,b,j,k,q,es;
@@ -1290,11 +1250,7 @@
subroutine to get augmented lagrange part of stiffness matrix
================================================================== */
-void get_aug_k(E,el,elt_k,level,m)
- struct All_variables *E;
- int el,m;
- double elt_k[24*24];
- int level;
+void get_aug_k(struct All_variables *E, int el, double elt_k[24*24], int level, int m)
{
int i,p[9],a,b,nodea,nodeb;
double Visc;
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_boundary_conditions.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_boundary_conditions.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -44,8 +44,7 @@
/* ========================================== */
-void full_velocity_boundary_conditions(E)
- struct All_variables *E;
+void full_velocity_boundary_conditions(struct All_variables *E)
{
void velocity_imp_vert_bc();
void velocity_apply_periodicapply_periodic_bcs();
@@ -131,10 +130,8 @@
/* ========================================== */
-void full_temperature_boundary_conditions(E)
- struct All_variables *E;
+void full_temperature_boundary_conditions(struct All_variables *E)
{
- void temperatures_conform_bcs();
void temperature_imposed_vert_bcs();
int j,lev,noz;
@@ -223,16 +220,14 @@
}
-static void velocity_apply_periodic_bcs(E)
- struct All_variables *E;
+static void velocity_apply_periodic_bcs(struct All_variables *E)
{
fprintf(E->fp,"Periodic boundary conditions\n");
return;
}
-static void temperature_apply_periodic_bcs(E)
- struct All_variables *E;
+static void temperature_apply_periodic_bcs(struct All_variables *E)
{
fprintf(E->fp,"Periodic temperature boundary conditions\n");
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_geometry_cartesian.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_geometry_cartesian.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_geometry_cartesian.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -31,8 +31,7 @@
#include "parsing.h"
-void full_set_2dc_defaults(E)
- struct All_variables *E;
+void full_set_2dc_defaults(struct All_variables *E)
{
E->mesh.nsd = 2;
@@ -41,8 +40,7 @@
}
-void full_set_2pt5dc_defaults(E)
- struct All_variables *E;
+void full_set_2pt5dc_defaults(struct All_variables *E)
{
E->mesh.nsd = 2;
@@ -50,8 +48,7 @@
}
-void full_set_3dc_defaults(E)
- struct All_variables *E;
+void full_set_3dc_defaults(struct All_variables *E)
{
E->mesh.nsd = 3;
@@ -59,8 +56,7 @@
}
-void full_set_3dsphere_defaults(E)
- struct All_variables *E;
+void full_set_3dsphere_defaults(struct All_variables *E)
{
E->mesh.nsd = 3;
E->mesh.dof = 3;
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_parallel_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_parallel_related.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_parallel_related.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -72,7 +72,7 @@
parallel_process_termination();
}
- E->sphere.caps_per_proc = max(1,E->sphere.caps*E->parallel.nprocz/E->parallel.nproc);
+ E->sphere.caps_per_proc = fmax(1,E->sphere.caps*E->parallel.nprocz/E->parallel.nproc);
if (E->sphere.caps_per_proc > 1) {
if (E->parallel.me==0) fprintf(stderr,"!!!! # caps per proc > 1 is not supported.\n \n");
@@ -326,8 +326,7 @@
exchange info across the boundaries
============================================ */
-void full_parallel_domain_boundary_nodes(E)
- struct All_variables *E;
+void full_parallel_domain_boundary_nodes(struct All_variables *E)
{
void parallel_process_termination();
@@ -502,8 +501,7 @@
static void face_eqn_node_to_pass(struct All_variables *, int, int, int, int);
static void line_eqn_node_to_pass(struct All_variables *, int, int, int, int, int, int);
-void full_parallel_communication_routs_v(E)
- struct All_variables *E;
+void full_parallel_communication_routs_v(struct All_variables *E)
{
int m,i,ii,j,k,l,node,el,elt,lnode,jj,doff,target;
@@ -759,8 +757,7 @@
assuming fault nodes are in the top row of processors
============================================ */
-void full_parallel_communication_routs_s(E)
- struct All_variables *E;
+void full_parallel_communication_routs_s(struct All_variables *E)
{
int i,ii,j,k,l,node,el,elt,lnode,jj,doff;
@@ -834,9 +831,7 @@
/* ================================================ */
/* ================================================ */
-static void face_eqn_node_to_pass(E,lev,m,npass,bd)
- struct All_variables *E;
- int lev,m,npass,bd;
+static void face_eqn_node_to_pass(struct All_variables *E, int lev, int m, int npass, int bd)
{
int jj,kk,node,doff;
const int dims=E->mesh.nsd;
@@ -859,9 +854,7 @@
/* ================================================ */
/* ================================================ */
-static void line_eqn_node_to_pass(E,lev,m,npass,num_node,offset,stride)
- struct All_variables *E;
- int lev,m,npass,num_node,offset,stride;
+static void line_eqn_node_to_pass(struct All_variables *E, int lev, int m, int npass, int num_node, int offset, int stride)
{
int jj,kk,node,doff;
const int dims=E->mesh.nsd;
@@ -902,10 +895,7 @@
by Tan2 7/21, 2003
================================================ */
-void full_exchange_id_d(E, U, lev)
- struct All_variables *E;
- double **U;
- int lev;
+void full_exchange_id_d(struct All_variables *E, double **U, int lev)
{
int ii,j,jj,m,k,kk,t_cap,idb,msginfo[8];
@@ -927,7 +917,7 @@
sizeofk = 0;
for (k=1;k<=E->parallel.TNUM_PASSz[lev];k++) {
kk = (1+E->parallel.NUM_NEQz[lev].pass[k])*sizeof(double);
- sizeofk = max(sizeofk, kk);
+ sizeofk = fmax(sizeofk, kk);
}
RV=(double *)malloc( sizeofk );
SV=(double *)malloc( sizeofk );
@@ -1018,10 +1008,7 @@
/* ================================================ */
/* ================================================ */
-static void exchange_node_d(E, U, lev)
- struct All_variables *E;
- double **U;
- int lev;
+static void exchange_node_d(struct All_variables *E, double **U, int lev)
{
int ii,j,jj,m,k,kk,t_cap,idb,msginfo[8];
@@ -1045,7 +1032,7 @@
idb= 0;
for (k=1;k<=E->parallel.TNUM_PASSz[lev];k++) {
sizeofk = (1+E->parallel.NUM_NODEz[lev].pass[k])*sizeof(double);
- idb = max(idb,sizeofk);
+ idb = fmax(idb,sizeofk);
}
RV=(double *)malloc( idb );
@@ -1142,10 +1129,7 @@
/* ================================================ */
/* ================================================ */
-static void exchange_node_f(E, U, lev)
- struct All_variables *E;
- float **U;
- int lev;
+static void exchange_node_f(struct All_variables *E, float **U, int lev)
{
int ii,j,jj,m,k,kk,t_cap,idb,msginfo[8];
@@ -1170,7 +1154,7 @@
idb= 0;
for (k=1;k<=E->parallel.TNUM_PASSz[lev];k++) {
sizeofk = (1+E->parallel.NUM_NODEz[lev].pass[k])*sizeof(float);
- idb = max(idb,sizeofk);
+ idb = fmax(idb,sizeofk);
}
RV=(float *)malloc( idb );
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_read_input_from_files.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_read_input_from_files.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -38,9 +38,7 @@
Open these files, read in results, and average if necessary
=========================================================================*/
-void full_read_input_files_for_timesteps(E,action,output)
- struct All_variables *E;
- int action, output;
+void full_read_input_files_for_timesteps(struct All_variables *E, int action, int output)
{
float find_age_in_MY();
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_sphere_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_sphere_related.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_sphere_related.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -168,11 +168,8 @@
double theta, fi, cost, sint, cosf, sinf, efac2,rfac;
double a, b;
double offset;
- double myatan();
- void even_divide_arc12();
-
- temp = max(E->mesh.noy, E->mesh.nox);
+ temp = fmax(E->mesh.noy, E->mesh.nox);
theta0 = (double *)malloc((temp+1)*sizeof(double));
fi0 = (double *)malloc((temp+1)*sizeof(double));
@@ -198,10 +195,10 @@
v4 = (double *)malloc((temp+1)*sizeof(double));
temp = E->mesh.noy * E->mesh.nox;
- px = malloc((temp+1)*sizeof(double));
- py = malloc((temp+1)*sizeof(double));
- qx = malloc((temp+1)*sizeof(double));
- qy = malloc((temp+1)*sizeof(double));
+ px = (double *)malloc((temp+1)*sizeof(double));
+ py = (double *)malloc((temp+1)*sizeof(double));
+ qx = (double *)malloc((temp+1)*sizeof(double));
+ qy = (double *)malloc((temp+1)*sizeof(double));
/* define the cap corners */
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_tracer_advection.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -160,22 +160,6 @@
/* Determine number of tracer quantities */
- /* advection_quantites - those needed for advection */
- E->trace.number_of_basic_quantities=12;
-
- /* extra_quantities - used for flavors, composition, etc. */
- /* (can be increased for additional science i.e. tracing chemistry */
-
- E->trace.number_of_extra_quantities = 0;
- if (E->trace.nflavors > 0)
- E->trace.number_of_extra_quantities += 1;
-
-
- E->trace.number_of_tracer_quantities =
- E->trace.number_of_basic_quantities +
- E->trace.number_of_extra_quantities;
-
-
/* Fixed positions in tracer array */
/* Flavor is always in extraq position 0 */
/* Current coordinates are always kept in basicq positions 0-5 */
@@ -184,22 +168,6 @@
/* Some error control regarding size of pointer arrays */
- if (E->trace.number_of_basic_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of basic in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of extraq in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_tracer_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
write_trace_instructions(E);
@@ -288,7 +256,6 @@
double *receive_z[13][3];
double *REC[13];
- void expand_tracer_arrays();
int icheck_that_processor_shell();
double CPU_time0();
@@ -312,7 +279,7 @@
fprintf(E->trace.fpt, "Entering lost_souls()\n");
- E->trace.istat_isend=E->trace.ilater[j];
+ E->trace.istat_isend=E->trace.later_tracer[j].size();
/** debug **
for (kk=1; kk<=E->trace.istat_isend; kk++) {
fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
@@ -327,13 +294,13 @@
/* initialize isend and ireceive */
/* # of neighbors in the horizontal plane */
- isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=E->trace.later_tracer[j].size()*(12+1);
for (kk=0;kk<=num_ngb;kk++) isend[j][kk]=0;
for (kk=0;kk<=num_ngb;kk++) ireceive[j][kk]=0;
/* Allocate Maximum Memory to Send Arrays */
- itemp_size=max(isize[j],1);
+ itemp_size=fmax(isize[j],1);
for (kk=0;kk<=num_ngb;kk++) {
if ((send[j][kk]=(double *)malloc(itemp_size*sizeof(double)))==NULL) {
@@ -407,9 +374,9 @@
/* Allocate memory in receive arrays */
for (ithatcap=0;ithatcap<=num_ngb;ithatcap++) {
- isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ isize[j]=ireceive[j][ithatcap]*(12+1);
- itemp_size=max(1,isize[j]);
+ itemp_size=fmax(1,isize[j]);
if ((receive[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
@@ -428,7 +395,7 @@
if (E->parallel.nprocz>1) {
ithatcap=ithiscap;
- isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ isize[j]=isend[j][ithatcap]*(12+1);
for (mm=0;mm<isize[j];mm++) {
receive[j][ithatcap][mm]=send[j][ithatcap][mm];
}
@@ -440,12 +407,12 @@
for (kk=1;kk<=num_ngb;kk++) {
idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
- isize[j]=isend[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=isend[j][kk]*(12+1);
MPI_Isend(send[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
- isize[j]=ireceive[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=ireceive[j][kk]*(12+1);
MPI_Irecv(receive[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
@@ -475,8 +442,8 @@
/* Allocate Memory for REC array */
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
+ isize[j]=isum[j]*(12+1);
+ isize[j]=fmax(isize[j],1);
if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
fflush(E->trace.fpt);
@@ -495,9 +462,9 @@
for (pp=0;pp<ireceive[j][ithatcap];pp++) {
irec[j]++;
- ipos=pp*E->trace.number_of_tracer_quantities;
+ ipos=pp*(12+1);
- for (mm=0;mm<E->trace.number_of_tracer_quantities;mm++) {
+ for (mm=0;mm<12+1;mm++) {
ipos2=ipos+mm;
REC[j][irec_position]=receive[j][ithatcap][ipos2];
@@ -520,8 +487,8 @@
/* (No dynamic reallocation of send_z necessary) */
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
- isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
+ isize[j]=itracers_subject_to_vertical_transport[j]*(12+1);
+ isize[j]=fmax(isize[j],1);
if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
@@ -546,7 +513,7 @@
num_tracers=irec[j];
for (kk=1;kk<=num_tracers;kk++) {
- ireceive_position=it*E->trace.number_of_tracer_quantities;
+ ireceive_position=it*(12+1);
it++;
irad=ireceive_position+2;
@@ -561,12 +528,12 @@
if (ival==1) {
- isend_position=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+ isend_position=isend_z[j][ivertical_neighbor]*(12+1);
isend_z[j][ivertical_neighbor]++;
- ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+ ilast_receiver_position=(irec[j]-1)*(12+1);
- for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++) {
+ for (mm=0;mm<12+1;mm++) {
ipos=ireceive_position+mm;
ipos2=isend_position+mm;
@@ -626,8 +593,8 @@
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
- isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
- isize[j]=max(isize[j],1);
+ isize[j]=ireceive_z[j][kk]*(12+1);
+ isize[j]=fmax(isize[j],1);
if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (t590)\n");
@@ -643,12 +610,12 @@
idestination_proc = E->parallel.PROCESSORz[lev].pass[kk];
- isize_send=isend_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize_send=isend_z[j][kk]*(12+1);
MPI_Isend(send_z[j][kk],isize_send,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
- isize_receive=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize_receive=ireceive_z[j][kk]*(12+1);
MPI_Irecv(receive_z[j][kk],isize_receive,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
@@ -670,7 +637,7 @@
isum[j]=isum[j]+irec[j];
- isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=isum[j]*(12+1);
if (isize[j]>0) {
if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL) {
@@ -686,11 +653,11 @@
for (kk=0;kk<ireceive_z[j][ivertical_neighbor];kk++) {
- irec_position=irec[j]*E->trace.number_of_tracer_quantities;
+ irec_position=irec[j]*(12+1);
irec[j]++;
- ireceive_position=kk*E->trace.number_of_tracer_quantities;
+ ireceive_position=kk*(12+1);
- for (mm=0;mm<E->trace.number_of_tracer_quantities;mm++) {
+ for (mm=0;mm<12+1;mm++) {
REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
}
}
@@ -711,41 +678,34 @@
for (kk=0;kk<irec[j];kk++) {
- E->trace.ntracers[j]++;
+ Tracer new_tracer;
- if (E->trace.ntracers[j]>(E->trace.max_ntracers[j]-5)) expand_tracer_arrays(E,j);
+ ireceive_position=kk*(12+1);
+ new_tracer.theta = REC[j][ipos+0];
+ new_tracer.phi = REC[j][ipos+1];
+ new_tracer.rad = REC[j][ipos+2];
+ new_tracer.x = REC[j][ipos+3];
+ new_tracer.y = REC[j][ipos+4];
+ new_tracer.z = REC[j][ipos+5];
+ new_tracer.x0 = REC[j][ipos+6];
+ new_tracer.y0 = REC[j][ipos+7];
+ new_tracer.z0 = REC[j][ipos+8];
+ new_tracer.Vx = REC[j][ipos+9];
+ new_tracer.Vy = REC[j][ipos+10];
+ new_tracer.Vz = REC[j][ipos+11];
+ new_tracer.flavor = REC[j][ipos+12];
- ireceive_position=kk*E->trace.number_of_tracer_quantities;
+ iel=(E->trace.iget_element)(E,j,-99,new_tracer.x,new_tracer.y,new_tracer.z,new_tracer.theta,new_tracer.phi,new_tracer.rad);
- for (mm=0;mm<E->trace.number_of_basic_quantities;mm++) {
- ipos=ireceive_position+mm;
-
- E->trace.basicq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
- }
- for (mm=0;mm<E->trace.number_of_extra_quantities;mm++) {
- ipos=ireceive_position+E->trace.number_of_basic_quantities+mm;
-
- E->trace.extraq[j][mm][E->trace.ntracers[j]]=REC[j][ipos];
- }
-
- theta=E->trace.basicq[j][0][E->trace.ntracers[j]];
- phi=E->trace.basicq[j][1][E->trace.ntracers[j]];
- rad=E->trace.basicq[j][2][E->trace.ntracers[j]];
- x=E->trace.basicq[j][3][E->trace.ntracers[j]];
- y=E->trace.basicq[j][4][E->trace.ntracers[j]];
- z=E->trace.basicq[j][5][E->trace.ntracers[j]];
-
-
- iel=(E->trace.iget_element)(E,j,-99,x,y,z,theta,phi,rad);
-
if (iel<1) {
fprintf(E->trace.fpt,"Error(lost souls) - element not here?\n");
- fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",x,y,z,theta,phi,rad);
+ fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",new_tracer.x,new_tracer.y,new_tracer.z,new_tracer.theta,new_tracer.phi,new_tracer.rad);
fflush(E->trace.fpt);
exit(10);
}
- E->trace.ielement[j][E->trace.ntracers[j]]=iel;
+ new_tracer.ielement = iel;
+ E->trace.tracers[j].push_back(new_tracer);
}
if(E->control.verbose){
@@ -778,22 +738,21 @@
{
const int j = 1;
int kk, pp;
- int numtracers, ithatcap, icheck;
+ int ithatcap, icheck;
int isend_position, ipos;
int lev = E->mesh.levmax;
double theta, phi, rad;
double x, y, z;
+ TracerList::iterator tr;
/* transfer tracers from rlater to send */
- numtracers=E->trace.ilater[j];
+ for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
+ rad=tr->rad;
+ x=tr->x;
+ y=tr->y;
+ z=tr->z;
- for (kk=1;kk<=numtracers;kk++) {
- rad=E->trace.rlater[j][2][kk];
- x=E->trace.rlater[j][3][kk];
- y=E->trace.rlater[j][4][kk];
- z=E->trace.rlater[j][5][kk];
-
/* first check same cap if nprocz>1 */
if (E->parallel.nprocz>1) {
@@ -829,12 +788,21 @@
/* assign tracer to send */
- isend_position=(isend[j][ithatcap]-1)*E->trace.number_of_tracer_quantities;
+ isend_position=(isend[j][ithatcap]-1)*(12+1);
- for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++) {
- ipos=isend_position+pp;
- send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
- }
+ send[j][ithatcap][isend_position+0]=tr->theta;
+ send[j][ithatcap][isend_position+1]=tr->phi;
+ send[j][ithatcap][isend_position+2]=tr->rad;
+ send[j][ithatcap][isend_position+3]=tr->x;
+ send[j][ithatcap][isend_position+4]=tr->y;
+ send[j][ithatcap][isend_position+5]=tr->z;
+ send[j][ithatcap][isend_position+6]=tr->x0;
+ send[j][ithatcap][isend_position+7]=tr->y0;
+ send[j][ithatcap][isend_position+8]=tr->z0;
+ send[j][ithatcap][isend_position+9]=tr->Vx;
+ send[j][ithatcap][isend_position+10]=tr->Vy;
+ send[j][ithatcap][isend_position+11]=tr->Vz;
+ send[j][ithatcap][isend_position+12]=tr->flavor;
} /* end kk, assigning tracers */
@@ -1304,8 +1272,8 @@
double *tmin;
double *tmax;
- double *fmin;
- double *fmax;
+ double *fmin_array;
+ double *fmax_array;
void sphere_to_cart();
@@ -1353,18 +1321,18 @@
theta=E->sx[j][1][kk];
phi=E->sx[j][2][kk];
- thetamax=max(thetamax,theta);
- thetamin=min(thetamin,theta);
+ thetamax=fmax(thetamax,theta);
+ thetamin=fmin(thetamin,theta);
}
/* expand range slightly (should take care of poles) */
thetamax=thetamax+expansion;
- thetamax=min(thetamax,M_PI);
+ thetamax=fmin(thetamax,M_PI);
thetamin=thetamin-expansion;
- thetamin=max(thetamin,0.0);
+ thetamin=fmax(thetamin,0.0);
/* Convert input data from degrees to radians */
@@ -1470,13 +1438,13 @@
fflush(E->trace.fpt);
exit(10);
}
- if ((fmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
+ if ((fmin_array=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
fflush(E->trace.fpt);
exit(10);
}
- if ((fmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
+ if ((fmax_array=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
fflush(E->trace.fpt);
@@ -1498,10 +1466,10 @@
theta=E->sx[j][1][node];
phi=E->sx[j][2][node];
- theta_min=min(theta_min,theta);
- theta_max=max(theta_max,theta);
- phi_min=min(phi_min,phi);
- phi_max=max(phi_max,phi);
+ theta_min=fmin(theta_min,theta);
+ theta_max=fmax(theta_max,theta);
+ phi_min=fmin(phi_min,phi);
+ phi_max=fmax(phi_max,phi);
}
/* add half difference to phi and expansion to theta to be safe */
@@ -1509,8 +1477,8 @@
theta_max=theta_max+expansion;
theta_min=theta_min-expansion;
- theta_max=min(M_PI,theta_max);
- theta_min=max(0.0,theta_min);
+ theta_max=fmin(M_PI,theta_max);
+ theta_min=fmax(0.0,theta_min);
half_diff=0.5*(phi_max-phi_min);
phi_max=phi_max+half_diff;
@@ -1527,8 +1495,8 @@
tmin[kk]=theta_min;
tmax[kk]=theta_max;
- fmin[kk]=phi_min;
- fmax[kk]=phi_max;
+ fmin_array[kk]=phi_min;
+ fmax_array[kk]=phi_max;
}
/* end looking through elements */
@@ -1565,7 +1533,7 @@
for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
{
pp=mm/elz;
- if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
+ if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin_array[pp]) && (phi<=fmax_array[pp]) )
{
ival=icheck_element_column(E,j,mm,x,y,z,rad);
if (ival>0)
@@ -1605,7 +1573,7 @@
for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
{
pp=mm/elz;
- if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
+ if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin_array[pp]) && (phi<=fmax_array[pp]) )
{
ival=icheck_element_column(E,j,mm,x,y,z,rad);
if (ival>0)
@@ -1630,8 +1598,8 @@
free(tmin);
free(tmax);
- free(fmin);
- free(fmax);
+ free(fmin_array);
+ free(fmax_array);
} /* end j */
@@ -1970,12 +1938,14 @@
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
+ /*
fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
E->trace.number_of_basic_quantities);
fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
E->trace.number_of_extra_quantities);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
E->trace.number_of_tracer_quantities);
+ */
/* analytical test */
@@ -2319,8 +2289,6 @@
double tiny, eps;
double x,y,z;
- double myatan();
-
/* make vectors from node to node */
makevector(v12,rnode2[1],rnode2[2],rnode2[3],rnode1[1],rnode1[2],rnode1[3]);
@@ -2866,7 +2834,7 @@
double u, v, cosc, theta_f, phi_f, dphi, cosd;
double *cost, *sint, *cosf, *sinf;
- if ((E->gnomonic = malloc((E->lmesh.nsf+1)*sizeof(struct CITCOM_GNOMONIC)))
+ if ((E->gnomonic = (CITCOM_GNOMONIC*)malloc((E->lmesh.nsf+1)*sizeof(struct CITCOM_GNOMONIC)))
== NULL) {
fprintf(stderr,"Error(define uv)-not enough memory(a)\n");
exit(10);
@@ -3066,9 +3034,7 @@
/* This function (and the 2 following) are used to test advection of tracers by assigning */
/* a test function (in "analytical_test_function"). */
-void analytical_test(E)
- struct All_variables *E;
-
+void analytical_test(struct All_variables *E)
{
#if 0
int kk,pp;
@@ -3325,15 +3291,7 @@
/*************** ANALYTICAL RUNGE KUTTE ******************/
/* */
-void analytical_runge_kutte(E,nsteps,dt,x0_s,x0_c,xf_s,xf_c,vec)
- struct All_variables *E;
- int nsteps;
- double dt;
- double *x0_c;
- double *x0_s;
- double *xf_c;
- double *xf_s;
- double *vec;
+void analytical_runge_kutte(struct All_variables *E, int nsteps, double dt, double *x0_s, double *x0_c, double *xf_s, double *xf_c, double *vec)
{
@@ -3445,12 +3403,7 @@
/* vel_c[2] => velocity in y direction */
/* vel_c[3] => velocity in z direction */
-void analytical_test_function(E,theta,phi,rad,vel_s,vel_c)
- struct All_variables *E;
- double theta,phi,rad;
- double *vel_s;
- double *vel_c;
-
+void analytical_test_function(struct All_variables *E, double theta, double phi, double rad, double *vel_s, double *vel_c)
{
double sint,sinf,cost,cosf;
Modified: mc/3D/CitcomS/branches/eheien/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Full_version_dependent.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Full_version_dependent.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -48,7 +48,7 @@
int m, int icap)
{
int i,lev;
- double t[4], myatan();
+ double t[4];
for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
for (i=1;i<=E->lmesh.NNO[lev];i++) {
@@ -78,8 +78,7 @@
================================================= */
-void full_node_locations(E)
- struct All_variables *E;
+void full_node_locations(struct All_variables *E)
{
int i,j,k,ii,lev;
double ro,dr,*rr,*RR,fo,tg;
Modified: mc/3D/CitcomS/branches/eheien/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/General_matrix_functions.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/General_matrix_functions.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -45,12 +45,7 @@
Iterative solver also using multigrid ........
=========================================================== */
-int solve_del2_u(E,d0,F,acc,high_lev)
- struct All_variables *E;
- double **d0;
- double **F;
- double acc;
- int high_lev;
+int solve_del2_u(struct All_variables *E, double **d0, double **F, double acc, int high_lev)
{
void assemble_del2_u();
void e_assemble_del2_u();
@@ -149,12 +144,7 @@
recursive multigrid function ....
================================= */
-double multi_grid(E,d1,F,acc,hl)
- struct All_variables *E;
- double **d1;
- double **F;
- double acc;
- int hl; /* higher level of two */
+double multi_grid(struct All_variables *E, double **d1, double **F, double acc, int hl)
{
double residual,AudotAu;
void interp_vector();
@@ -304,13 +294,7 @@
#ifndef USE_CUDA
-double conj_grad(E,d0,F,acc,cycles,level)
- struct All_variables *E;
- double **d0;
- double **F;
- double acc;
- int *cycles;
- int level;
+double conj_grad(struct All_variables *E, double **d0, double **F, double acc, int *cycles, int level)
{
double *r0[NCS],*r1[NCS],*r2[NCS];
double *z0[NCS],*z1[NCS],*z2[NCS];
@@ -431,14 +415,7 @@
versions
=========================================================================================*/
-void element_gauss_seidel(E,d0,F,Ad,acc,cycles,level,guess)
- struct All_variables *E;
- double **d0;
- double **F,**Ad;
- double acc;
- int *cycles;
- int level;
- int guess;
+void element_gauss_seidel(struct All_variables *E, double **d0, double **F, double **Ad, double acc, int *cycles, int level, int guess)
{
int count,i,j,k,l,m,ns,nc,d,steps,loc;
int p1,p2,p3,q1,q2,q3;
@@ -603,14 +580,7 @@
time (Jacobi at a node). It does the job though.
============================================================================ */
-void gauss_seidel(E,d0,F,Ad,acc,cycles,level,guess)
- struct All_variables *E;
- double **d0;
- double **F,**Ad;
- double acc;
- int *cycles;
- int level;
- int guess;
+void gauss_seidel(struct All_variables *E, double **d0, double **F, double **Ad, double acc, int *cycles, int level, int guess)
{
int count,i,j,k,l,m,ns,steps;
@@ -758,9 +728,7 @@
/* Fast (conditional) determinant for 3x3 or 2x2 ... otherwise calls general routine */
-double determinant(A,n)
- double A[4][4];
- int n;
+double determinant(double A[4][4], int n)
{ double gen_determinant();
@@ -781,12 +749,9 @@
-double cofactor(A,i,j,n)
- double A[4][4];
- int i,j,n;
+double cofactor(double A[4][4], int i, int j, int n)
{ int k,l,p,q;
- double determinant();
double B[4][4]; /* because of recursive behaviour of det/cofac, need to use
new copy of B at each 'n' level of this routine */
Modified: mc/3D/CitcomS/branches/eheien/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Ggrd_handling.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Ggrd_handling.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -33,7 +33,7 @@
*/
#ifdef USE_GZDIR
#include <zlib.h>
-gzFile *gzdir_output_open(char *,char *);
+//gzFile *gzdir_output_open(char *,char *);
#endif
Modified: mc/3D/CitcomS/branches/eheien/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Global_operations.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Global_operations.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -44,10 +44,7 @@
aren't & another method is required.
=============================================== */
-void remove_horiz_ave(E,X,H,store_or_not)
- struct All_variables *E;
- double **X, *H;
- int store_or_not;
+void remove_horiz_ave(struct All_variables *E, double **X, double *H, int store_or_not)
{
int m,i,j,k,n,nox,noz,noy;
@@ -83,9 +80,7 @@
}
-void return_horiz_ave(E,X,H)
- struct All_variables *E;
- double **X, *H;
+void return_horiz_ave(struct All_variables *E, double **X, double *H)
{
const int dims = E->mesh.nsd;
int m,i,j,k,d,nint,noz,nox,noy,el,elz,elx,ely,j1,j2,i1,i2,k1,k2,nproc;
@@ -163,9 +158,7 @@
return;
}
-void return_horiz_ave_f(E,X,H)
- struct All_variables *E;
- float **X, *H;
+void return_horiz_ave_f(struct All_variables *E, float **X, float *H)
{
const int dims = E->mesh.nsd;
int m,i,j,k,d,nint,noz,nox,noy,el,elz,elx,ely,j1,j2,i1,i2,k1,k2,nproc;
@@ -250,9 +243,7 @@
/* however here, elemental horizontal averages are given rather than */
/* nodal averages. Also note, here is average per element */
-void return_elementwise_horiz_ave(E,X,H)
- struct All_variables *E;
- double **X, *H;
+void return_elementwise_horiz_ave(struct All_variables *E, double **X, double *H)
{
int m,i,j,k,d,noz,noy,el,elz,elx,ely,nproc;
@@ -312,11 +303,7 @@
return;
}
-float return_bulk_value(E,Z,average)
- struct All_variables *E;
- float **Z;
- int average;
-
+float return_bulk_value(struct All_variables *E, float **Z, int average)
{
int n,i,j,k,el,m;
float volume,integral,volume1,integral1;
@@ -356,11 +343,7 @@
/* when integer average =0, integral is returned. */
-double return_bulk_value_d(E,Z,average)
- struct All_variables *E;
- double **Z;
- int average;
-
+double return_bulk_value_d(struct All_variables *E, double **Z, int average)
{
int n,i,j,el,m;
double volume,integral,volume1,integral1;
@@ -394,9 +377,7 @@
}
/* ================================================== */
-float find_max_horizontal(E,Tmax)
-struct All_variables *E;
-float Tmax;
+float find_max_horizontal(struct All_variables *E, float Tmax)
{
float ttmax;
@@ -406,10 +387,7 @@
}
/* ================================================== */
-void sum_across_surface(E,data,total)
-struct All_variables *E;
-float *data;
-int total;
+void sum_across_surface(struct All_variables *E, float *data, int total)
{
int j,d;
float *temp;
@@ -430,9 +408,7 @@
/* ================================================== */
/* ================================================== */
-void sum_across_surf_sph1(E,sphc,sphs)
-struct All_variables *E;
-float *sphc,*sphs;
+void sum_across_surf_sph1(struct All_variables *E, float *sphc, float *sphs)
{
int jumpp,total,j,d;
float *sphcs,*temp;
@@ -466,11 +442,7 @@
/* ================================================== */
-float global_fvdot(E,A,B,lev)
- struct All_variables *E;
- float **A,**B;
- int lev;
-
+float global_fvdot(struct All_variables *E, float **A, float **B, int lev)
{
int m,i,neq;
float prod, temp,temp1;
@@ -499,11 +471,7 @@
}
-double kineticE_radial(E,A,lev)
- struct All_variables *E;
- double **A;
- int lev;
-
+double kineticE_radial(struct All_variables *E, double **A, int lev)
{
int m,i,neq;
double prod, temp,temp1;
@@ -531,11 +499,7 @@
return (prod);
}
-double global_vdot(E,A,B,lev)
- struct All_variables *E;
- double **A,**B;
- int lev;
-
+double global_vdot(struct All_variables *E, double **A, double **B, int lev)
{
int m,i,neq;
double prod, temp,temp1;
@@ -562,11 +526,7 @@
}
-double global_pdot(E,A,B,lev)
- struct All_variables *E;
- double **A,**B;
- int lev;
-
+double global_pdot(struct All_variables *E, double **A, double **B, int lev)
{
int i,m,npno;
double prod, temp;
@@ -656,11 +616,7 @@
}
-double global_tdot_d(E,A,B,lev)
- struct All_variables *E;
- double **A,**B;
- int lev;
-
+double global_tdot_d(struct All_variables *E, double **A, double **B, int lev)
{
int i,nno,m;
double prod, temp;
@@ -681,11 +637,7 @@
return (prod);
}
-float global_tdot(E,A,B,lev)
- struct All_variables *E;
- float **A,**B;
- int lev;
-
+float global_tdot(struct All_variables *E, float **A, float **B, int lev)
{
int i,nno,m;
float prod, temp;
@@ -706,18 +658,14 @@
}
-float global_fmin(E,a)
- struct All_variables *E;
- float a;
+float global_fmin(struct All_variables *E, double a)
{
float temp;
MPI_Allreduce(&a, &temp,1,MPI_FLOAT,MPI_MIN,E->parallel.world);
return (temp);
}
-double global_dmax(E,a)
- struct All_variables *E;
- double a;
+double global_dmax(struct All_variables *E, double a)
{
double temp;
MPI_Allreduce(&a, &temp,1,MPI_DOUBLE,MPI_MAX,E->parallel.world);
@@ -725,18 +673,14 @@
}
-float global_fmax(E,a)
- struct All_variables *E;
- float a;
+float global_fmax(struct All_variables *E, float a)
{
float temp;
MPI_Allreduce(&a, &temp,1,MPI_FLOAT,MPI_MAX,E->parallel.world);
return (temp);
}
-double Tmaxd(E,T)
- struct All_variables *E;
- double **T;
+double Tmaxd(struct All_variables *E, double **T)
{
double global_dmax(),temp,temp1;
int i,m;
@@ -744,16 +688,14 @@
temp = -10.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
- temp = max(T[m][i],temp);
+ temp = fmax(T[m][i],temp);
temp1 = global_dmax(E,temp);
return (temp1);
}
-float Tmax(E,T)
- struct All_variables *E;
- float **T;
+float Tmax(struct All_variables *E, float **T)
{
float global_fmax(),temp,temp1;
int i,m;
@@ -761,17 +703,14 @@
temp = -10.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++)
- temp = max(T[m][i],temp);
+ temp = fmax(T[m][i],temp);
temp1 = global_fmax(E,temp);
return (temp1);
}
-double vnorm_nonnewt(E,dU,U,lev)
- struct All_variables *E;
- double **dU,**U;
- int lev;
+double vnorm_nonnewt(struct All_variables *E, double **dU, double **U, int lev)
{
double temp1,temp2,dtemp,temp;
int a,e,i,m,node;
@@ -803,9 +742,7 @@
}
-void sum_across_depth_sph1(E,sphc,sphs)
- struct All_variables *E;
- float *sphc,*sphs;
+void sum_across_depth_sph1(struct All_variables *E, float *sphc, float *sphs)
{
int jumpp,total,j;
@@ -888,7 +825,6 @@
void remove_rigid_rot(struct All_variables *E)
{
void velo_from_element_d();
- double myatan();
double wx, wy, wz, v_theta, v_phi, cos_t,sin_t,sin_f, cos_f,frd;
double vx[9], vy[9], vz[9];
double r, t, f, efac,tg;
Modified: mc/3D/CitcomS/branches/eheien/lib/Initial_temperature.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Initial_temperature.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Initial_temperature.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -288,7 +288,7 @@
}
/* Truncate the temperature to be within (0,1). */
/* This might not be desirable in some situations. */
- E->T[m][i] = max(0.0,min(g,1.0));
+ E->T[m][i] = fmax(0.0,fmin(g,1.0));
}
}
fclose (fp);
Modified: mc/3D/CitcomS/branches/eheien/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Instructions.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Instructions.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -623,7 +623,7 @@
heat flux files? */
if(E->output.write_q_files){ /* make sure those get written at
least as often as velocities */
- E->output.write_q_files = min(E->output.write_q_files,E->control.record_every);
+ E->output.write_q_files = fmin(E->output.write_q_files,E->control.record_every);
}
@@ -898,9 +898,7 @@
common to all problems follow ...
=================================== */
-void allocate_common_vars(E)
- struct All_variables *E;
-
+void allocate_common_vars(struct All_variables *E)
{
void set_up_nonmg_aliases();
int m,n,snel,nsf,elx,ely,nox,noy,noz,nno,nel,npno,lim;
@@ -955,8 +953,8 @@
/* nodal mass */
E->NMass[j] = (double *) malloc((nno+1)*sizeof(double));
- nxyz = max(nox*noz,nox*noy);
- nxyz = 2*max(nxyz,noz*noy);
+ nxyz = fmax(nox*noz,nox*noy);
+ nxyz = 2*fmax(nxyz,noz*noy);
E->sien[j] = (struct SIEN *) malloc((nxyz+2)*sizeof(struct SIEN));
E->surf_element[j] = (int *) malloc((nxyz+2)*sizeof(int));
@@ -1051,9 +1049,9 @@
E->VI[i][j] = (float *) malloc((nno+1)*sizeof(float));
E->NODE[i][j] = (unsigned int *)malloc((nno+1)*sizeof(unsigned int));
- nxyz = max(nox*noz,nox*noy);
- nxyz = 2*max(nxyz,noz*noy);
- nozl = max(noy,nox*2);
+ nxyz = fmax(nox*noz,nox*noy);
+ nxyz = 2*fmax(nxyz,noz*noy);
+ nozl = fmax(noy,nox*2);
@@ -1134,9 +1132,7 @@
/* ========================================================= */
-void allocate_velocity_vars(E)
- struct All_variables *E;
-
+void allocate_velocity_vars(struct All_variables *E)
{
int m,n,i,j,k,l;
@@ -1196,8 +1192,7 @@
/* ========================================================= */
-void global_default_values(E)
- struct All_variables *E;
+void global_default_values(struct All_variables *E)
{
/* FIRST: values which are not changed routinely by the user */
@@ -1296,8 +1291,7 @@
/* =============================================================
============================================================= */
-void check_bc_consistency(E)
- struct All_variables *E;
+void check_bc_consistency(struct All_variables *E)
{ int i,j,lev;
@@ -1336,9 +1330,7 @@
}
-void set_up_nonmg_aliases(E,j)
- struct All_variables *E;
- int j;
+void set_up_nonmg_aliases(struct All_variables *E, int j)
{ /* Aliases for functions only interested in the highest mg level */
@@ -1363,9 +1355,7 @@
return; }
-void report(E,string)
- struct All_variables *E;
- char * string;
+void report(struct All_variables *E, char * string)
{ if(E->control.verbose && E->parallel.me==0)
{ fprintf(stderr,"%s\n",string);
fflush(stderr);
@@ -1373,9 +1363,7 @@
return;
}
-void record(E,string)
- struct All_variables *E;
- char * string;
+void record(struct All_variables *E, char * string)
{ if(E->control.verbose && E->fp)
{ fprintf(E->fp,"%s\n",string);
fflush(E->fp);
@@ -1395,8 +1383,7 @@
/* This function is replaced by CitcomS.Components.IC.launch()*/
-void common_initial_fields(E)
- struct All_variables *E;
+void common_initial_fields(struct All_variables *E)
{
void initial_pressure();
void initial_velocity();
@@ -1412,8 +1399,7 @@
/* ========================================== */
-void initial_pressure(E)
- struct All_variables *E;
+void initial_pressure(struct All_variables *E)
{
int i,m;
report(E,"Initialize pressure field");
@@ -1425,8 +1411,7 @@
return;
}
-void initial_velocity(E)
- struct All_variables *E;
+void initial_velocity(struct All_variables *E)
{
int i,m;
report(E,"Initialize velocity field");
Modified: mc/3D/CitcomS/branches/eheien/lib/Interuption.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Interuption.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Interuption.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -29,10 +29,9 @@
#include <stdio.h>
#include <stdlib.h>
+#include "global_defs.h"
#include "interuption.h"
-void parallel_process_termination();
-
int Emergency_stop;
Modified: mc/3D/CitcomS/branches/eheien/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Makefile.am 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Makefile.am 2011-07-14 00:06:30 UTC (rev 18754)
@@ -60,6 +60,7 @@
# static library
libCitcomS_a_CFLAGS = $(AM_CFLAGS) # hack for automake
+libCitcomS_a_CXXFLAGS = $(AM_CXXFLAGS) # hack for automake
libCitcomS_a_SOURCES = $(sources)
# shared library (libtool)
@@ -129,7 +130,7 @@
Stokes_flow_Incomp.c \
Topo_gravity.c \
tracer_defs.h \
- Tracer_setup.c \
+ Tracer_setup.cpp \
viscosity_descriptions.h \
Viscosity_structures.c \
Full_boundary_conditions.c \
Modified: mc/3D/CitcomS/branches/eheien/lib/Mineral_physics_models.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Mineral_physics_models.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Mineral_physics_models.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -176,10 +176,10 @@
/* reference model (PREM) */
- rhor = malloc((E->lmesh.noz+1) * sizeof(double));
- vpr = malloc((E->lmesh.noz+1) * sizeof(double));
- vsr = malloc((E->lmesh.noz+1) * sizeof(double));
- depthkm = malloc((E->lmesh.noz+1) * sizeof(double));
+ rhor = (double *)malloc((E->lmesh.noz+1) * sizeof(double));
+ vpr = (double *)malloc((E->lmesh.noz+1) * sizeof(double));
+ vsr = (double *)malloc((E->lmesh.noz+1) * sizeof(double));
+ depthkm = (double *)malloc((E->lmesh.noz+1) * sizeof(double));
for(nz=1; nz<=E->lmesh.noz; nz++) {
get_prem(E->sx[m][3][nz], &vpr[nz], &vsr[nz], &rhor[nz]);
Modified: mc/3D/CitcomS/branches/eheien/lib/Nodal_mesh.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Nodal_mesh.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Nodal_mesh.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -36,8 +36,7 @@
/* get nodal spherical velocities from the solution vector */
-void v_from_vector(E)
- struct All_variables *E;
+void v_from_vector(struct All_variables *E)
{
int m,node;
const int nno = E->lmesh.nno;
@@ -59,8 +58,7 @@
return;
}
-void assign_v_to_vector(E)
- struct All_variables *E;
+void assign_v_to_vector(struct All_variables *E)
{
int m,node;
const int nno = E->lmesh.nno;
@@ -75,8 +73,7 @@
return;
}
-void v_from_vector_pseudo_surf(E)
- struct All_variables *E;
+void v_from_vector_pseudo_surf(struct All_variables *E)
{
int m,node;
@@ -118,10 +115,7 @@
return;
}
/* cartesian velocities within element, single prec version */
-void velo_from_element(E,VV,m,el,sphere_key)
- struct All_variables *E;
- float VV[4][9];
- int el,m,sphere_key;
+void velo_from_element(struct All_variables *E, float VV[4][9], int m, int el, int sphere_key)
{
int a, node;
@@ -159,10 +153,7 @@
}
/* double prec version */
-void velo_from_element_d(E,VV,m,el,sphere_key)
- struct All_variables *E;
- double VV[4][9];
- int el,m,sphere_key;
+void velo_from_element_d(struct All_variables *E, double VV[4][9], int m, int el, int sphere_key)
{
int a, node;
@@ -202,11 +193,7 @@
}
-void p_to_nodes(E,P,PN,lev)
- struct All_variables *E;
- double **P;
- float **PN;
- int lev;
+void p_to_nodes(struct All_variables *E, double **P, float **PN, int lev)
{ int e,element,node,j,m;
@@ -236,10 +223,7 @@
interpolate the viscosity from element integration points to nodes
*/
-void visc_from_gint_to_nodes(E,VE,VN,lev)
- struct All_variables *E;
- float **VE,**VN;
- int lev;
+void visc_from_gint_to_nodes(struct All_variables *E, float **VE, float **VN, int lev)
{
int m,e,i,j,k,n,off,lim;
const int nsd=E->mesh.nsd;
@@ -275,10 +259,7 @@
interpolate viscosity from nodes to element integration points
*/
-void visc_from_nodes_to_gint(E,VN,VE,lev)
- struct All_variables *E;
- float **VE,**VN;
- int lev;
+void visc_from_nodes_to_gint(struct All_variables *E, float **VN, float **VE, int lev)
{
int m,e,i,j,k,n,off;
@@ -309,10 +290,7 @@
visc_from_gint_to_ele(E,E->EVI[lv],viscU,lv)
*/
-void visc_from_gint_to_ele(E,VE,VN,lev)
- struct All_variables *E;
- float **VE,**VN;
- int lev;
+void visc_from_gint_to_ele(struct All_variables *E, float **VE, float **VN, int lev)
{
int m,e,i,j,k,n,off;
const int nsd=E->mesh.nsd;
@@ -341,10 +319,7 @@
*/
-void visc_from_ele_to_gint(E,VN,VE,lev)
- struct All_variables *E;
- float **VE,**VN;
- int lev;
+void visc_from_ele_to_gint(struct All_variables *E, float **VN, float **VE, int lev)
{
int m,e,i,j,k,n,off;
const int nsd=E->mesh.nsd;
Modified: mc/3D/CitcomS/branches/eheien/lib/Output.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Output.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Output.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -562,9 +562,9 @@
double *rho, *vp, *vs;
const int len = E->lmesh.nno;
- rho = malloc(len * sizeof(double));
- vp = malloc(len * sizeof(double));
- vs = malloc(len * sizeof(double));
+ rho = (double *)malloc(len * sizeof(double));
+ vp = (double *)malloc(len * sizeof(double));
+ vs = (double *)malloc(len * sizeof(double));
if(rho==NULL || vp==NULL || vs==NULL) {
fprintf(stderr, "Error while allocating memory\n");
abort();
@@ -659,28 +659,29 @@
int i, j, n, ncolumns;
char output_file[255];
FILE *fp1;
+ TracerArray::iterator ta;
+ TracerList::iterator tr;
+
sprintf(output_file,"%s.tracer.%d.%d", E->control.data_file,
E->parallel.me, cycles);
fp1 = output_open(output_file, "w");
- ncolumns = 3 + E->trace.number_of_extra_quantities;
+ //ncolumns = 3 + E->trace.number_of_extra_quantities;
+ ncolumns = 3 + 1;
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- fprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
+ for(ta=E->trace.tracers.begin();ta!=E->trace.tracers.end();++ta) {
+ fprintf(fp1,"%d %d %d %.5e\n", cycles, ta->size(),
ncolumns, E->monitor.elapsed_time);
- for(n=1;n<=E->trace.ntracers[j];n++) {
+ for(tr=ta->begin();tr!=ta->end();++tr) {
/* write basic quantities (coordinate) */
- fprintf(fp1,"%.12e %.12e %.12e",
- E->trace.basicq[j][0][n],
- E->trace.basicq[j][1][n],
- E->trace.basicq[j][2][n]);
+ fprintf(fp1,"%.12e %.12e %.12e %.12e",
+ tr->theta,
+ tr->phi,
+ tr->rad,
+ tr->flavor);
- /* write extra quantities */
- for (i=0; i<E->trace.number_of_extra_quantities; i++) {
- fprintf(fp1," %.12e", E->trace.extraq[j][i][n]);
- }
fprintf(fp1, "\n");
}
Modified: mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Output_gzdir.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -1078,29 +1078,29 @@
int i, j, n, ncolumns;
char output_file[255];
gzFile *fp1;
+ TracerArray::iterator ta;
+ TracerList::iterator tr;
snprintf(output_file,255,"%s/%d/tracer.%d.%d.gz",
E->control.data_dir,cycles,
E->parallel.me, cycles);
fp1 = gzdir_output_open(output_file,"w");
- ncolumns = 3 + E->trace.number_of_extra_quantities;
+ ncolumns = 3 + 1;
+ //ncolumns = 3 + E->trace.number_of_extra_quantities;
- for(j=1;j<=E->sphere.caps_per_proc;j++) {
- gzprintf(fp1,"%d %d %d %.5e\n", cycles, E->trace.ntracers[j],
+ for(ta=E->trace.tracers.begin();ta!=E->trace.tracers.end();++ta) {
+ gzprintf(fp1,"%d %d %d %.5e\n", cycles, ta->size(),
ncolumns, E->monitor.elapsed_time);
- for(n=1;n<=E->trace.ntracers[j];n++) {
+ for(tr=ta->begin();tr!=ta->end();++tr) {
/* write basic quantities (coordinate) */
- gzprintf(fp1,"%9.5e %9.5e %9.5e",
- E->trace.basicq[j][0][n],
- E->trace.basicq[j][1][n],
- E->trace.basicq[j][2][n]);
+ gzprintf(fp1,"%9.5e %9.5e %9.5e %9.5e",
+ tr->theta,
+ tr->phi,
+ tr->rad,
+ tr->flavor);
- /* write extra quantities */
- for (i=0; i<E->trace.number_of_extra_quantities; i++) {
- gzprintf(fp1," %9.5e", E->trace.extraq[j][i][n]);
- }
gzprintf(fp1, "\n");
}
@@ -1514,8 +1514,8 @@
/* this should not be called with (i,i,size i) */
void be_flipit(void *d, void *s, size_t len)
{
- unsigned char *dest = d;
- unsigned char *src = s;
+ unsigned char *dest = (unsigned char *)d;
+ unsigned char *src = (unsigned char *)s;
src += len - 1;
for (; len; len--)
*dest++ = *src--;
Modified: mc/3D/CitcomS/branches/eheien/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Pan_problem_misc_functions.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Pan_problem_misc_functions.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -80,9 +80,7 @@
}
-void unique_copy_file(E,name,comment)
- struct All_variables *E;
- char *name, *comment;
+void unique_copy_file(struct All_variables *E, char *name, char *comment)
{
char unique_name[500];
char command[600];
@@ -337,8 +335,7 @@
my version of arc tan
=================================================*/
-double myatan(y,x)
- double y,x;
+double myatan(double y, double x)
{
double fi;
Modified: mc/3D/CitcomS/branches/eheien/lib/Parsing.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Parsing.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Parsing.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -84,12 +84,9 @@
int interpret_control_string();
-void setup_parser(E,filename)
- struct All_variables *E;
- char *filename;
+void setup_parser(struct All_variables *E, char *filename)
{
void unique_copy_file();
- void add_to_parameter_list();
FILE * fp;
char *pl,*pn,*pv;
@@ -175,9 +172,7 @@
}
-void shutdown_parser(E)
- struct All_variables *E;
-
+void shutdown_parser(struct All_variables *E)
{
if(ARGLIST != NULL) free(ARGLIST);
if(ARGBUF != NULL) free(ARGBUF);
@@ -188,13 +183,11 @@
/* add an entry to arglist, expanding memory */
-void add_to_parameter_list(name,value)
- char *name, *value; /* if necessary */
+void add_to_parameter_list(char *name, char *value)
{
struct arglist *alptr;
int len;
char *ptr;
- int compute_parameter_hash_table();
/* check arglist memory */
if(NLIST >= LISTMAX)
@@ -233,8 +226,7 @@
NLIST++;
}
-int compute_parameter_hash_table(s)
- char *s;
+int compute_parameter_hash_table(char *s)
{ int h;
h= s[0];
@@ -251,12 +243,7 @@
return(h);
}
-int input_int(name,value,interpret,m)
- char *name;
- int *value;
- char *interpret;
- int m;
-
+int input_int(char *name, int *value, char *interpret, int m)
{
struct arglist *alptr;
int h, found;
@@ -320,11 +307,8 @@
return(found);
}
-int input_string(name,value,Default,m) /* in the case of a string default=NULL forces input */
- char *name;
- char *value;
- char *Default;
- int m;
+/* in the case of a string default=NULL forces input */
+int input_string(char *name, char *value, char *Default, int m)
{
char *sptr;
struct arglist *alptr;
@@ -374,11 +358,8 @@
return(found);
}
-int input_boolean(name,value,interpret,m) /* supports name=on/off too */
- char *name;
- int *value;
- char *interpret;
- int m;
+/* supports name=on/off too */
+int input_boolean(char *name, int *value, char *interpret, int m)
{
char *sptr;
struct arglist *alptr;
@@ -443,11 +424,7 @@
return(found);
}
-int input_float(name,value,interpret,m)
- char *name;
- float *value;
- char *interpret;
- int m;
+int input_float(char *name, float *value, char *interpret, int m)
{ char *sptr;
struct arglist *alptr;
@@ -511,11 +488,7 @@
return(found);
}
-int input_double(name,value,interpret,m)
- char *name;
- double *value;
- char *interpret;
- int m;
+int input_double(char *name, double *value, char *interpret, int m)
{ char *sptr;
struct arglist *alptr;
@@ -630,11 +603,7 @@
-int input_char_vector(name,number,value,m)
- char *name;
- int number;
- char *value; /* comma-separated list of ints */
- int m;
+int input_char_vector(char *name, int number, char *value, int m)
{ char *sptr;
struct arglist *alptr;
@@ -679,11 +648,7 @@
return(found);
}
-int input_float_vector(name,number,value,m)
- char *name;
- int number;
- float *value; /* comma-separated list of floats */
- int m;
+int input_float_vector(char *name, int number, float *value, int m)
{ char *sptr;
struct arglist *alptr;
@@ -731,11 +696,7 @@
return(found);
}
-int input_double_vector(name,number,value,m)
- char *name;
- int number;
- double *value; /* comma-separated list of floats */
- int m;
+int input_double_vector(char *name, int number, double *value, int m)
{ char *sptr;
struct arglist *alptr;
@@ -786,10 +747,7 @@
The function strtok does not work on linux machine
*/
-int interpret_control_string(interpret,essential,Default,minvalue,maxvalue)
- char *interpret;
- int *essential;
- double *Default,*minvalue,*maxvalue;
+int interpret_control_string(char *interpret, int *essential, double *Default, double *minvalue, double *maxvalue)
{ char *substring;
Modified: mc/3D/CitcomS/branches/eheien/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Problem_related.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Problem_related.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -38,8 +38,7 @@
read velocity vectors at the top surface from files
=========================================================================*/
-void read_velocity_boundary_from_file(E)
- struct All_variables *E;
+void read_velocity_boundary_from_file(struct All_variables *E)
{
(E->solver.read_input_files_for_timesteps)(E,1,1); /* read velocity(1) and output(1) */
return;
@@ -51,8 +50,7 @@
laterally varying rayleigh number in the top layers
*/
-void read_rayleigh_from_file(E)
- struct All_variables *E;
+void read_rayleigh_from_file(struct All_variables *E)
{
(E->solver.read_input_files_for_timesteps)(E,4,1); /* read Rayleigh number for top layers */
return;
@@ -62,8 +60,7 @@
construct material array
=========================================================================*/
-void read_mat_from_file(E)
- struct All_variables *E;
+void read_mat_from_file(struct All_variables *E)
{
(E->solver.read_input_files_for_timesteps)(E,3,1); /* read element material(3) and output(1) */
return;
@@ -73,8 +70,7 @@
read temperature at the top surface from files
=========================================================================*/
-void read_temperature_boundary_from_file(E)
- struct All_variables *E;
+void read_temperature_boundary_from_file(struct All_variables *E)
{
(E->solver.read_input_files_for_timesteps)(E,5,1); /* read temperature(5) and output(1) */
return;
@@ -85,8 +81,7 @@
Open restart file to get initial elapsed time, or calculate the right value
=========================================================================*/
-void get_initial_elapsed_time(E)
- struct All_variables *E;
+void get_initial_elapsed_time(struct All_variables *E)
{
FILE *fp;
int ll, mm,rezip;
@@ -136,8 +131,7 @@
Sets the elapsed time to zero, if desired.
=========================================================================*/
-void set_elapsed_time(E)
- struct All_variables *E;
+void set_elapsed_time(struct All_variables *E)
{
if (E->control.zero_elapsed_time) /* set elapsed_time to zero */
@@ -151,8 +145,7 @@
run, if desired.
=========================================================================*/
-void set_starting_age(E)
- struct All_variables *E;
+void set_starting_age(struct All_variables *E)
{
/* remember start_age is in MY */
if (E->control.reset_startage)
@@ -168,9 +161,7 @@
making ages SMALLER!
=========================================================================*/
- float find_age_in_MY(E)
-
- struct All_variables *E;
+float find_age_in_MY(struct All_variables *E)
{
float age_in_MY, e_4;
Modified: mc/3D/CitcomS/branches/eheien/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Process_buoyancy.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Process_buoyancy.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -50,8 +50,7 @@
Surface heat flux
=================== */
-void heat_flux(E)
- struct All_variables *E;
+void heat_flux(struct All_variables *E)
{
int m,e,el,i,j,node,lnode,nz;
float *flux[NCS],*SU[NCS],*RU[NCS];
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_boundary_conditions.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_boundary_conditions.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_boundary_conditions.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -39,15 +39,14 @@
void assign_internal_bc(struct All_variables * ,int);
static void velocity_apply_periodic_bcs();
static void temperature_apply_periodic_bcs();
-static void velocity_refl_vert_bc();
-static void temperature_refl_vert_bc();
+static void velocity_refl_vert_bc(struct All_variables *E);
+static void temperature_refl_vert_bc(struct All_variables *E);
void read_temperature_boundary_from_file(struct All_variables *);
void read_velocity_boundary_from_file(struct All_variables *);
/* ========================================== */
-void regional_velocity_boundary_conditions(E)
- struct All_variables *E;
+void regional_velocity_boundary_conditions(struct All_variables *E)
{
void velocity_imp_vert_bc();
void renew_top_velocity_boundary();
@@ -144,12 +143,10 @@
/* ========================================== */
-void regional_temperature_boundary_conditions(E)
- struct All_variables *E;
+void regional_temperature_boundary_conditions(struct All_variables *E)
{
void temperature_imposed_vert_bcs();
void temperature_lith_adj();
- void temperatures_conform_bcs();
int j,lev,noz;
lev = E->mesh.levmax;
@@ -196,8 +193,7 @@
/* ========================================== */
-static void velocity_refl_vert_bc(E)
- struct All_variables *E;
+static void velocity_refl_vert_bc(struct All_variables *E)
{
int m,i,j,ii,jj;
int node1,node2;
@@ -339,8 +335,7 @@
return;
}
-static void temperature_refl_vert_bc(E)
- struct All_variables *E;
+static void temperature_refl_vert_bc(struct All_variables *E)
{
int i,j,m;
int node1,node2;
@@ -436,8 +431,7 @@
}
-static void velocity_apply_periodic_bcs(E)
- struct All_variables *E;
+static void velocity_apply_periodic_bcs(struct All_variables *E)
{
int n1,n2,level;
int i,j,ii,jj;
@@ -448,8 +442,7 @@
return;
}
-static void temperature_apply_periodic_bcs(E)
- struct All_variables *E;
+static void temperature_apply_periodic_bcs(struct All_variables *E)
{
const int dims=E->mesh.nsd;
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_geometry_cartesian.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_geometry_cartesian.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_geometry_cartesian.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -29,8 +29,7 @@
#include "parsing.h"
-void regional_set_2dc_defaults(E)
- struct All_variables *E;
+void regional_set_2dc_defaults(struct All_variables *E)
{
E->mesh.nsd = 2;
@@ -39,8 +38,7 @@
}
-void regional_set_2pt5dc_defaults(E)
- struct All_variables *E;
+void regional_set_2pt5dc_defaults(struct All_variables *E)
{
E->mesh.nsd = 2;
@@ -48,8 +46,7 @@
}
-void regional_set_3dc_defaults(E)
- struct All_variables *E;
+void regional_set_3dc_defaults(struct All_variables *E)
{
E->mesh.nsd = 3;
@@ -57,8 +54,7 @@
}
-void regional_set_3dsphere_defaults(E)
- struct All_variables *E;
+void regional_set_3dsphere_defaults(struct All_variables *E)
{
E->mesh.nsd = 3;
E->mesh.dof = 3;
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_parallel_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_parallel_related.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_parallel_related.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -294,8 +294,7 @@
exchange info across the boundaries
============================================ */
-void regional_parallel_domain_boundary_nodes(E)
- struct All_variables *E;
+void regional_parallel_domain_boundary_nodes(struct All_variables *E)
{
void parallel_process_termination();
@@ -443,8 +442,7 @@
assuming fault nodes are in the top row of processors
============================================ */
-void regional_parallel_communication_routs_v(E)
- struct All_variables *E;
+void regional_parallel_communication_routs_v(struct All_variables *E)
{
int m,i,ii,j,k,l,node,el,elt,lnode,jj,doff,target_cap;
@@ -619,8 +617,7 @@
assuming fault nodes are in the top row of processors
============================================ */
-void regional_parallel_communication_routs_s(E)
- struct All_variables *E;
+void regional_parallel_communication_routs_s(struct All_variables *E)
{
int i,ii,j,k,l,node,el,elt,lnode,jj,doff;
@@ -741,10 +738,7 @@
by Tan2 7/21, 2003
================================================ */
-void regional_exchange_id_d(E, U, lev)
- struct All_variables *E;
- double **U;
- int lev;
+void regional_exchange_id_d(struct All_variables *E, double **U, int lev)
{
int ii,j,jj,m,k;
@@ -791,10 +785,7 @@
/* ================================================ */
/* ================================================ */
-static void exchange_node_d(E, U, lev)
- struct All_variables *E;
- double **U;
- int lev;
+static void exchange_node_d(struct All_variables *E, double **U, int lev)
{
int ii,j,jj,m,k;
@@ -840,10 +831,7 @@
/* ================================================ */
/* ================================================ */
-static void exchange_node_f(E, U, lev)
- struct All_variables *E;
- float **U;
- int lev;
+static void exchange_node_f(struct All_variables *E, float **U, int lev)
{
int ii,j,jj,m,k;
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_read_input_from_files.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_read_input_from_files.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -38,9 +38,7 @@
Open these files, read in results, and average if necessary
=========================================================================*/
-void regional_read_input_files_for_timesteps(E,action,output)
- struct All_variables *E;
- int action, output;
+void regional_read_input_files_for_timesteps(struct All_variables *E, int action, int output)
{
float find_age_in_MY();
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_sphere_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_sphere_related.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_sphere_related.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -32,9 +32,7 @@
#include "element_definitions.h"
#include "global_defs.h"
-void regional_coord_of_cap(E,m,icap)
- struct All_variables *E;
- int icap,m;
+void regional_coord_of_cap(struct All_variables *E, int m, int icap)
{
int i,j,k,lev,temp,elx,ely,nox,noy,noz,node,nodes;
@@ -69,7 +67,7 @@
nprocyl=E->parallel.nprocy;
nproczl=E->parallel.nprocz;
nnproc=nprocyl*nprocxl*nproczl;
- temp = max(E->mesh.NOY[E->mesh.levmax],E->mesh.NOX[E->mesh.levmax]);
+ temp = fmax(E->mesh.NOY[E->mesh.levmax],E->mesh.NOX[E->mesh.levmax]);
/* define the cap corners */
E->sphere.cap[1].theta[1] = E->control.theta_min;
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_tracer_advection.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -49,7 +49,7 @@
static void make_mesh_ijk(struct All_variables *E);
static void put_lost_tracers(struct All_variables *E,
int *send_size, double *send,
- int kk, int j);
+ Tracer send_tracer, int j);
static void put_found_tracers(struct All_variables *E,
int recv_size, double *recv,
int j);
@@ -92,48 +92,12 @@
/* This parameter specifies how close a tracer can get to the boundary */
E->trace.box_cushion=0.00001;
- /* Determine number of tracer quantities */
-
- /* advection_quantites - those needed for advection */
- E->trace.number_of_basic_quantities=12;
-
- /* extra_quantities - used for flavors, composition, etc. */
- /* (can be increased for additional science i.e. tracing chemistry */
-
- E->trace.number_of_extra_quantities = 0;
- if (E->trace.nflavors > 0)
- E->trace.number_of_extra_quantities += 1;
-
-
- E->trace.number_of_tracer_quantities =
- E->trace.number_of_basic_quantities +
- E->trace.number_of_extra_quantities;
-
-
/* Fixed positions in tracer array */
/* Flavor is always in extraq position 0 */
/* Current coordinates are always kept in basicq positions 0-5 */
/* Other positions may be used depending on science being done */
- /* Some error control regarding size of pointer arrays */
-
- if (E->trace.number_of_basic_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of basic in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_extra_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of extraq in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
- if (E->trace.number_of_tracer_quantities>99) {
- fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in tracer_defs.h\n");
- fflush(E->trace.fpt);
- parallel_process_termination();
- }
-
write_trace_instructions(E);
/* The bounding box of neiboring processors */
@@ -209,12 +173,14 @@
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
+ /*
fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
E->trace.number_of_basic_quantities);
fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
E->trace.number_of_extra_quantities);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
E->trace.number_of_tracer_quantities);
+ */
@@ -634,7 +600,6 @@
double *send[2];
double *recv[2];
- void expand_tracer_arrays();
int icheck_that_processor_shell();
int ipass;
@@ -645,7 +610,7 @@
double CPU_time0();
double begin_time = CPU_time0();
- E->trace.istat_isend = E->trace.ilater[j];
+ E->trace.istat_isend = E->trace.later_tracer[j].size();
/* the bounding box */
for (d=0; d<E->mesh.nsd; d++) {
@@ -689,8 +654,8 @@
/* Allocate Maximum Memory to Send Arrays */
- max_send_size = max(2*E->trace.ilater[j], E->trace.ntracers[j]/100);
- itemp_size = max_send_size * E->trace.number_of_tracer_quantities;
+ max_send_size = fmax(2*E->trace.later_tracer[j].size(), E->trace.tracers[j].size()/100);
+ itemp_size = max_send_size * (12+1);
if ((send[0] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
@@ -706,70 +671,74 @@
}
+ TracerList::iterator tr;
+
for (d=0; d<E->mesh.nsd; d++) {
- int original_size = E->trace.ilater[j];
+ int original_size = E->trace.later_tracer[j].size();
int idb;
- int kk = 1;
int isend[2], irecv[2];
isend[0] = isend[1] = 0;
-
/* move out-of-bound tracers to send array */
- while (kk<=E->trace.ilater[j]) {
+ for (tr=E->trace.later_tracer[j].begin();tr!=E->trace.later_tracer[j].end();++tr) {
double coord;
/* Is the tracer within the bounds in the d-th dimension */
- coord = E->trace.rlater[j][d][kk];
+ switch(d) {
+ case 0:
+ coord=tr->theta;
+ break;
+ case 1:
+ coord=tr->phi;
+ break;
+ case 2:
+ coord=tr->rad;
+ break;
+ }
if (coord < bounds[d][0]) {
- put_lost_tracers(E, &(isend[0]), send[0], kk, j);
+ put_lost_tracers(E, &(isend[0]), send[0], *tr, j);
+ tr = E->trace.later_tracer[j].erase(tr);
}
else if (coord >= bounds[d][1]) {
- put_lost_tracers(E, &(isend[1]), send[1], kk, j);
+ put_lost_tracers(E, &(isend[1]), send[1], *tr, j);
+ tr = E->trace.later_tracer[j].erase(tr);
}
else {
/* check next tracer */
- kk++;
}
/* reallocate send if size too small */
- if ((isend[0] > max_send_size - 5) ||
- (isend[1] > max_send_size - 5)) {
+ if ((isend[0] > max_send_size - 5) || (isend[1] > max_send_size - 5)) {
isize = max_send_size + max_send_size/4 + 10;
- itemp_size = isize * E->trace.number_of_tracer_quantities;
+ itemp_size = isize * (12+1);
- if ((send[0] = (double *)realloc(send[0],
- itemp_size*sizeof(double)))
- == NULL) {
+ if ((send[0] = (double *)realloc(send[0], itemp_size*sizeof(double))) == NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (s4)\n");
fflush(E->trace.fpt);
exit(10);
}
- if ((send[1] = (double *)realloc(send[1],
- itemp_size*sizeof(double)))
- == NULL) {
+ if ((send[1] = (double *)realloc(send[1], itemp_size*sizeof(double))) == NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (s5)\n");
fflush(E->trace.fpt);
exit(10);
}
- fprintf(E->trace.fpt,"Expanding physical memory of send to "
- "%d from %d\n",
+ fprintf(E->trace.fpt,"Expanding physical memory of send to %d from %d\n",
isize, max_send_size);
max_send_size = isize;
}
+ }
- } /* end of while kk */
-
/* check the total # of tracers is conserved */
- if ((isend[0] + isend[1] + E->trace.ilater[j]) != original_size) {
+ if ((isend[0] + isend[1] + E->trace.later_tracer[j].size()) != original_size) {
fprintf(E->trace.fpt, "original_size: %d, rlater_size: %d, "
"send_size: %d\n",
- original_size, E->trace.ilater[j], kk);
+ original_size, (int)E->trace.later_tracer[j].size(), kk);
}
@@ -827,8 +796,8 @@
/* Allocate memory in receive arrays */
for (i=0; i<2; i++) {
- isize = irecv[i] * E->trace.number_of_tracer_quantities;
- itemp_size = max(1, isize);
+ isize = irecv[i] * (12+1);
+ itemp_size = fmax(1, isize);
if ((recv[i] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
@@ -846,11 +815,11 @@
kk = d*2 + i + 1;
target_rank = ngbr_rank[kk];
if (target_rank >= 0) {
- isize = isend[i] * E->trace.number_of_tracer_quantities;
+ isize = isend[i] * (12+1);
MPI_Isend(send[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
- isize = irecv[i] * E->trace.number_of_tracer_quantities;
+ isize = irecv[i] * (12+1);
MPI_Irecv(recv[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
@@ -889,14 +858,16 @@
/* rlater should be empty by now */
- if (E->trace.ilater[j] > 0) {
+ if (E->trace.later_tracer[j].size() > 0) {
fprintf(E->trace.fpt, "Error(regional_lost_souls) lost tracers\n");
+ /*
for (kk=1; kk<=E->trace.ilater[j]; kk++) {
fprintf(E->trace.fpt, "lost #%d xx=(%e, %e, %e)\n", kk,
E->trace.rlater[j][0][kk],
E->trace.rlater[j][1][kk],
E->trace.rlater[j][2][kk]);
}
+ */
fflush(E->trace.fpt);
exit(10);
}
@@ -914,27 +885,30 @@
static void put_lost_tracers(struct All_variables *E,
int *send_size, double *send,
- int kk, int j)
+ Tracer send_tracer, int j)
{
int ilast_tracer, isend_position, ipos;
int pp;
/* move the tracer from rlater to send */
- isend_position = (*send_size) * E->trace.number_of_tracer_quantities;
+ isend_position = (*send_size) * (12+1);
- for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
- ipos = isend_position + pp;
- send[ipos] = E->trace.rlater[j][pp][kk];
- }
+ send[isend_position+0] = send_tracer.theta;
+ send[isend_position+1] = send_tracer.phi;
+ send[isend_position+2] = send_tracer.rad;
+ send[isend_position+3] = send_tracer.x;
+ send[isend_position+4] = send_tracer.y;
+ send[isend_position+5] = send_tracer.z;
+ send[isend_position+6] = send_tracer.x0;
+ send[isend_position+7] = send_tracer.y0;
+ send[isend_position+8] = send_tracer.z0;
+ send[isend_position+9] = send_tracer.Vx;
+ send[isend_position+10] = send_tracer.Vy;
+ send[isend_position+11] = send_tracer.Vz;
+ send[isend_position+12] = send_tracer.flavor;
+
(*send_size)++;
- /* eject the tracer from rlater */
- ilast_tracer = E->trace.ilater[j];
- for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
- E->trace.rlater[j][pp][kk] = E->trace.rlater[j][pp][ilast_tracer];
- }
- E->trace.ilater[j]--;
-
return;
}
@@ -947,25 +921,35 @@
int recv_size, double *recv,
int j)
{
- void expand_tracer_arrays();
- void expand_later_array();
int icheck_processor_shell();
int kk, pp;
- int ipos, ilast, inside, iel;
- double theta, phi, rad;
+ int ipos, inside, iel;
for (kk=0; kk<recv_size; kk++) {
- ipos = kk * E->trace.number_of_tracer_quantities;
- theta = recv[ipos];
- phi = recv[ipos + 1];
- rad = recv[ipos + 2];
+ Tracer new_tracer;
+ ipos = kk * (12+1);
+
+ /* found the element */
+ new_tracer.theta = recv[ipos+0];
+ new_tracer.phi = recv[ipos+1];
+ new_tracer.rad = recv[ipos+2];
+ new_tracer.x = recv[ipos+3];
+ new_tracer.y = recv[ipos+4];
+ new_tracer.z = recv[ipos+5];
+ new_tracer.x0 = recv[ipos+6];
+ new_tracer.y0 = recv[ipos+7];
+ new_tracer.z0 = recv[ipos+8];
+ new_tracer.Vx = recv[ipos+9];
+ new_tracer.Vy = recv[ipos+10];
+ new_tracer.Vz = recv[ipos+11];
+
/* check whether this tracer is inside this proc */
/* check radius first, since it is cheaper */
- inside = icheck_processor_shell(E, j, rad);
+ inside = icheck_processor_shell(E, j, new_tracer.rad);
if (inside == 1)
- inside = regional_icheck_cap(E, 0, theta, phi, rad, rad);
+ inside = regional_icheck_cap(E, 0, new_tracer.theta, new_tracer.phi, new_tracer.rad, new_tracer.rad);
else
inside = 0;
@@ -977,58 +961,21 @@
/**/
if (inside) {
+ iel = regional_iget_element(E, j, -99, 0, 0, 0, new_tracer.theta, new_tracer.phi, new_tracer.rad);
- E->trace.ntracers[j]++;
- ilast = E->trace.ntracers[j];
-
- if (E->trace.ntracers[j] > (E->trace.max_ntracers[j]-5))
- expand_tracer_arrays(E, j);
-
- for (pp=0; pp<E->trace.number_of_basic_quantities; pp++)
- E->trace.basicq[j][pp][ilast] = recv[ipos+pp];
-
- ipos += E->trace.number_of_basic_quantities;
- for (pp=0; pp<E->trace.number_of_extra_quantities; pp++)
- E->trace.extraq[j][pp][ilast] = recv[ipos+pp];
-
-
- /* found the element */
- iel = regional_iget_element(E, j, -99, 0, 0, 0, theta, phi, rad);
-
if (iel<1) {
- fprintf(E->trace.fpt, "Error(regional lost souls) - "
- "element not here?\n");
- fprintf(E->trace.fpt, "theta, phi, rad: %f %f %f\n",
- theta, phi, rad);
+ fprintf(E->trace.fpt, "Error(regional lost souls) - element not here?\n");
+ fprintf(E->trace.fpt, "theta, phi, rad: %f %f %f\n", new_tracer.theta, new_tracer.phi, new_tracer.rad);
fflush(E->trace.fpt);
exit(10);
}
- E->trace.ielement[j][ilast] = iel;
+ new_tracer.ielement = iel;
+ E->trace.tracers[j].push_back(new_tracer);
}
else {
- if (E->trace.ilatersize[j]==0) {
-
- E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
-
- for (kk=0;kk<E->trace.number_of_tracer_quantities;kk++) {
- if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL) {
- fprintf(E->trace.fpt,"AKM(put_found_tracers)-no memory (%d)\n",kk);
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- } /* end first particle initiating memory allocation */
-
- E->trace.ilater[j]++;
- ilast = E->trace.ilater[j];
-
- if (E->trace.ilater[j] > (E->trace.ilatersize[j]-5))
- expand_later_array(E, j);
-
- for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++)
- E->trace.rlater[j][pp][ilast] = recv[ipos+pp];
+ E->trace.later_tracer[j].push_back(new_tracer);
} /* end of if-else */
/** debug **
Modified: mc/3D/CitcomS/branches/eheien/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Regional_version_dependent.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Regional_version_dependent.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -43,8 +43,7 @@
================================================= */
-void regional_node_locations(E)
- struct All_variables *E;
+void regional_node_locations(struct All_variables *E)
{
int i,j,k,lev;
double ro,dr,*rr,*RR,fo;
Modified: mc/3D/CitcomS/branches/eheien/lib/Shape_functions.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Shape_functions.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Shape_functions.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -41,10 +41,8 @@
Function creating shape_fn data in form of a structure
=======================================================*/
-void construct_shape_functions(E)
- struct All_variables *E;
+void construct_shape_functions(struct All_variables *E)
{
- double lpoly(),lpolydash();
int i,j,k,d,dd;
int remapj,remapk;
@@ -172,9 +170,10 @@
return; }
-double lpoly(p,y)
- int p; /* selects lagrange polynomial , 1d: node p */
- double y; /* coordinate in given direction to evaluate poly */
+//int p; /* selects lagrange polynomial , 1d: node p */
+//double y; /* coordinate in given direction to evaluate poly */
+
+double lpoly(int p, double y)
{
double value;
@@ -193,9 +192,7 @@
return(value);
}
-double lpolydash(p,y)
- int p;
- double y;
+double lpolydash(int p, double y)
{
double value;
switch (p)
Modified: mc/3D/CitcomS/branches/eheien/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Size_does_matter.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Size_does_matter.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -34,8 +34,7 @@
double theta_g(double , struct All_variables *);
#endif
-void twiddle_thumbs(yawn)
- struct All_variables *yawn;
+void twiddle_thumbs(struct All_variables *yawn)
// int scratch_groin;
{ /* Do nothing, just sit back and relax.
@@ -51,8 +50,6 @@
static void form_rtf_bc(int k, double x[4],
double rtf[4][9], double bc[4][4])
{
- double myatan();
-
rtf[3][k] = 1.0/sqrt(x[1]*x[1]+x[2]*x[2]+x[3]*x[3]); /* 1/r */
rtf[1][k] = acos(x[3]*rtf[3][k]); /* theta */
rtf[2][k] = myatan(x[2],x[1]); /* phi */
@@ -76,9 +73,6 @@
{
int i,j,k,d,e;
double jacobian;
- double determinant();
- double cofactor(),myatan();
- void form_rtf_bc();
struct Shape_function_dx LGNx;
@@ -209,8 +203,6 @@
int i, k, d;
double x[4];
- double myatan();
-
const int dims = E->mesh.nsd;
const int ends = ENODES3D;
const int vpts = VPOINTS3D;
@@ -239,8 +231,6 @@
int i, k, d;
double x[4];
- double myatan();
-
const int dims = E->mesh.nsd;
const int ends = ENODES3D;
const int ppts = PPOINTS3D;
@@ -306,14 +296,12 @@
/* ======================================================================
====================================================================== */
-void construct_surf_det (E)
- struct All_variables *E;
+void construct_surf_det (struct All_variables *E)
{
int m,i,k,d,e,es,el;
double jacobian;
- double determinant();
double cofactor();
const int oned = onedvpoints[E->mesh.nsd];
@@ -363,7 +351,6 @@
int m,i,k,d,e,es,el,side;
double jacobian;
- double determinant();
double cofactor();
const int oned = onedvpoints[E->mesh.nsd];
@@ -418,16 +405,11 @@
/* ======================================================================
====================================================================== */
-void get_global_1d_shape_fn(E,el,GM,dGammax,top,m)
- struct All_variables *E;
- int el,top,m;
- struct Shape_function1 *GM;
- struct Shape_function1_dA *dGammax;
+void get_global_1d_shape_fn(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dA *dGammax, int top, int m)
{
int ii,i,k,d,e;
double jacobian;
- double determinant();
const int oned = onedvpoints[E->mesh.nsd];
@@ -457,16 +439,11 @@
/* ======================================================================
====================================================================== */
-void get_global_1d_shape_fn_L(E,el,GM,dGammax,top,m)
- struct All_variables *E;
- int el,top,m;
- struct Shape_function1 *GM;
- struct Shape_function1_dA *dGammax;
+void get_global_1d_shape_fn_L(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dA *dGammax, int top, int m)
{
int ii,i,k,d,e,node;
double jacobian;
- double determinant();
const int oned = onedvpoints[E->mesh.nsd];
@@ -527,17 +504,11 @@
/* ======================================================================
For calculating pressure boundary term --- Choi, 11/13/02
====================================================================== */
-void get_global_side_1d_shape_fn(E,el,GM,GMx,dGamma,side,m)
- struct All_variables *E;
- int el,side,m;
- struct Shape_function1 *GM;
- struct Shape_function1_dx *GMx;
- struct Shape_function_side_dA *dGamma;
+void get_global_side_1d_shape_fn(struct All_variables *E, int el, struct Shape_function1 *GM, struct Shape_function1_dx *GMx, struct Shape_function_side_dA *dGamma, int side, int m)
{
int i,k,d,e;
double jacobian;
- double determinant();
const int oned = onedvpoints[E->mesh.nsd];
double xx[4][5],dxda[4][4];
@@ -571,7 +542,7 @@
struct CCX *ccx,int lev,int m,int pressure)
{
int a,i,j,k,d,lnode;
- double cofactor(),myatan();
+ double cofactor();
double x[4],u[4][4],ux[3][4][4],ua[4][4];
double costt,cosff,sintt,sinff,rr,tt,ff;
@@ -713,7 +684,7 @@
int lev,int m,int pressure,int side)
{
int a,aa,i,j,k,d,lnode;
- double cofactor(),myatan();
+ double cofactor();
double x[4],u[4][4],ux[3][4][4],ua[4][4];
double costt,cosff,sintt,sinff,rr,tt,ff;
@@ -850,11 +821,10 @@
/* ======================================= */
-void construct_c3x3matrix(E)
- struct All_variables *E;
+void construct_c3x3matrix(struct All_variables *E)
{
int m,a,i,j,k,d,es,el,nel_surface,lev,lnode;
- double cofactor(),myatan();
+ double cofactor();
double x[4],u[4][4],ux[3][4][4],ua[4][4];
double costt,cosff,sintt,sinff,rr,tt,ff;
@@ -1009,7 +979,7 @@
{
int m,node,i,nint,e,lev;
int n[9], nz;
- double myatan(),area,centre[4],temp[9],temp2[9],dx1,dx2,dx3;
+ double area,centre[4],temp[9],temp2[9],dx1,dx2,dx3;
const int vpts=vpoints[E->mesh.nsd];
@@ -1046,7 +1016,7 @@
E->ECO[lev][m][e].centre[3] = dx3;
/* delta(theta) of this element */
- dx1 = max( fabs(E->SX[lev][m][1][n[3]]-E->SX[lev][m][1][n[1]]),
+ dx1 = fmax( fabs(E->SX[lev][m][1][n[3]]-E->SX[lev][m][1][n[1]]),
fabs(E->SX[lev][m][1][n[2]]-E->SX[lev][m][1][n[4]]) );
/* length of this element in the theta-direction */
@@ -1055,15 +1025,15 @@
/* delta(phi) of this element */
dx1 = fabs(E->SX[lev][m][2][n[3]]-E->SX[lev][m][2][n[1]]);
if (dx1>M_PI)
- dx1 = min(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) + 2.0*M_PI -
- max(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) ;
+ dx1 = fmin(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) + 2.0*M_PI -
+ fmax(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) ;
dx2 = fabs(E->SX[lev][m][2][n[2]]-E->SX[lev][m][2][n[4]]);
if (dx2>M_PI)
- dx2 = min(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) + 2.0*M_PI -
- max(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) ;
+ dx2 = fmin(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) + 2.0*M_PI -
+ fmax(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) ;
- dx2 = max(dx1,dx2);
+ dx2 = fmax(dx1,dx2);
/* length of this element in the phi-direction */
E->ECO[lev][m][e].size[2] = dx2*E->ECO[lev][m][e].centre[3]
Modified: mc/3D/CitcomS/branches/eheien/lib/Solver_conj_grad.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Solver_conj_grad.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Solver_conj_grad.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -28,11 +28,8 @@
#include "element_definitions.h"
#include "global_defs.h"
-void set_cg_defaults(E)
- struct All_variables *E;
-{ void assemble_forces_iterative();
- void solve_constrained_flow_iterative();
- void cg_allocate_vars();
+void set_cg_defaults(struct All_variables *E)
+{
E->control.CONJ_GRAD = 1;
E->build_forcing_term = assemble_forces_iterative;
@@ -43,8 +40,7 @@
return;
}
-void cg_allocate_vars(E)
- struct All_variables *E;
+void cg_allocate_vars(struct All_variables *E)
{
/* Nothing required ONLY by conj-grad stuff */
/* printf("here here\n"); */
@@ -53,8 +49,7 @@
}
-void assemble_forces_iterative(E)
- struct All_variables *E;
+void assemble_forces_iterative(struct All_variables *E)
{
int i;
Modified: mc/3D/CitcomS/branches/eheien/lib/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Solver_multigrid.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Solver_multigrid.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -32,12 +32,8 @@
#include "anisotropic_viscosity.h"
#endif
-void set_mg_defaults(E)
- struct All_variables *E;
-{ void assemble_forces_iterative();
- void solve_constrained_flow_iterative();
- void mg_allocate_vars();
-
+void set_mg_defaults(struct All_variables *E)
+{
E->control.NMULTIGRID = 1;
E->build_forcing_term = assemble_forces_iterative;
E->solve_stokes_problem = solve_constrained_flow_iterative;
@@ -47,8 +43,7 @@
return;
}
-void mg_allocate_vars(E)
- struct All_variables *E;
+void mg_allocate_vars(struct All_variables *E)
{
return;
@@ -60,11 +55,9 @@
just dropping values at shared grid points.
===================================================== */
-void inject_scalar(E,start_lev,AU,AD)
- struct All_variables *E;
- int start_lev;
- float **AU,**AD; /* data on upper/lower mesh */
+// float **AU,**AD; /* data on upper/lower mesh */
+void inject_scalar(struct All_variables *E, int start_lev, float **AU, float **AD)
{
int i,m,el,node_coarse,node_fine,sl_minus,eqn,eqn_coarse;
void gather_to_1layer_node ();
@@ -91,11 +84,9 @@
return;
}
-void inject_vector(E,start_lev,AU,AD)
- struct All_variables *E;
- int start_lev;
- double **AU,**AD; /* data on upper/lower mesh */
+// double **AU,**AD; /* data on upper/lower mesh */
+void inject_vector(struct All_variables *E, int start_lev, double **AU, double **AD)
{
int i,j,m,el,node_coarse,node_fine,sl_minus,eqn_fine,eqn_coarse;
@@ -132,11 +123,7 @@
just dropping values at shared grid points.
===================================================== */
-void un_inject_vector(E,start_lev,AD,AU)
-
- struct All_variables *E;
- int start_lev;
- double **AU,**AD; /* data on upper/lower mesh */
+void un_inject_vector(struct All_variables *E, int start_lev, double **AD, double **AU)
{
int i,m;
int el,node,node_plus;
@@ -184,11 +171,7 @@
======================================================================================= */
-void interp_vector(E,start_lev,AD,AU)
-
- struct All_variables *E;
- int start_lev;
- double **AD,**AU; /* data on upper/lower mesh */
+void interp_vector(struct All_variables *E, int start_lev, double **AD, double **AU)
{
void un_inject_vector();
void fill_in_gaps();
@@ -229,9 +212,7 @@
levels in the problem. (no gaps for vbcs)
============================================== */
-void project_viscosity(E)
- struct All_variables *E;
-
+void project_viscosity(struct All_variables *E)
{
int lv,i,j,k,m,sl_minus;
void inject_scalar();
@@ -350,11 +331,7 @@
}
/* ==================================================== */
-void inject_scalar_e(E,start_lev,AU,AD)
-
- struct All_variables *E;
- int start_lev;
- float **AU,**AD; /* data on upper/lower mesh */
+void inject_scalar_e(struct All_variables *E, int start_lev, float **AU, float **AD)
{
int i,j,m;
int el,node,e;
@@ -384,11 +361,7 @@
}
/* ==================================================== */
-void project_scalar_e(E,start_lev,AU,AD)
-
- struct All_variables *E;
- int start_lev;
- float **AU,**AD; /* data on upper/lower mesh */
+void project_scalar_e(struct All_variables *E, int start_lev, float **AU, float **AD)
{
int i,j,m;
int el,node,e;
@@ -423,11 +396,7 @@
}
/* ==================================================== */
-void project_scalar(E,start_lev,AU,AD)
-
- struct All_variables *E;
- int start_lev;
- float **AU,**AD; /* data on upper/lower mesh */
+void project_scalar(struct All_variables *E, int start_lev, float **AU, float **AD)
{
int i,j,m;
int el,node,node1;
@@ -476,11 +445,7 @@
/* this is prefered scheme with averages */
-void project_vector(E,start_lev,AU,AD,ic)
-
- struct All_variables *E;
- int start_lev,ic;
- double **AU,**AD; /* data on upper/lower mesh */
+void project_vector(struct All_variables *E, int start_lev, double **AU, double **AD, int ic)
{
int i,j,m;
int el,node1,node,e1;
@@ -553,10 +518,7 @@
}
/* ================================================= */
- void from_xyz_to_rtf(E,level,xyz,rtf)
- struct All_variables *E;
- int level;
- double **rtf,**xyz;
+ void from_xyz_to_rtf(struct All_variables *E, int level, double **xyz, double **rtf)
{
int i,j,m,eqn1,eqn2,eqn3;
@@ -585,10 +547,7 @@
}
/* ================================================= */
- void from_rtf_to_xyz(E,level,rtf,xyz)
- struct All_variables *E;
- int level;
- double **rtf,**xyz;
+ void from_rtf_to_xyz(struct All_variables *E, int level, double **rtf, double **xyz)
{
int i,j,m,eqn1,eqn2,eqn3;
@@ -618,10 +577,7 @@
}
/* ========================================================== */
- void fill_in_gaps(E,temp,level)
- struct All_variables *E;
- int level;
- double **temp;
+ void fill_in_gaps(struct All_variables *E, double **temp, int level)
{
int i,j,k,m;
Modified: mc/3D/CitcomS/branches/eheien/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Sphere_harmonics.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Sphere_harmonics.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -12,9 +12,7 @@
/* ======================================================================
====================================================================== */
-void set_sphere_harmonics(E)
- struct All_variables *E;
-
+void set_sphere_harmonics(struct All_variables *E)
{
int m,node,ll,mm,i,j;
@@ -100,9 +98,7 @@
/* =========================================================
expand the field TG into spherical harmonics
========================================================= */
-void sphere_expansion(E,TG,sphc,sphs)
- struct All_variables *E;
- float **TG,*sphc,*sphs;
+void sphere_expansion(struct All_variables *E, float **TG, float *sphc, float *sphs)
{
int el,nint,d,p,i,m,j,es,mm,ll,rand();
void sum_across_surf_sph1();
@@ -193,11 +189,8 @@
/* ==================================================*/
/* ==================================================*/
-static void compute_sphereh_table(E)
- struct All_variables *E;
+static void compute_sphereh_table(struct All_variables *E)
{
- double modified_plgndr_a();
-
int m,node,ll,mm,i,j,p;
double t,f,mmf;
Modified: mc/3D/CitcomS/branches/eheien/lib/Sphere_util.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Sphere_util.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Sphere_util.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -37,11 +37,9 @@
anglewise evenly.
=================================================*/
-void even_divide_arc12(elx,x1,y1,z1,x2,y2,z2,theta,fi)
- double x1,y1,z1,x2,y2,z2,*theta,*fi;
- int elx;
+void even_divide_arc12(int elx, double x1, double y1, double z1, double x2, double y2, double z2, double *theta, double *fi)
{
- double dx,dy,dz,xx,yy,zz,myatan();
+ double dx,dy,dz,xx,yy,zz;
int j, nox;
nox = elx+1;
@@ -64,13 +62,11 @@
compute angle and area
================================================*/
-void compute_angle_surf_area (E)
- struct All_variables *E;
+void compute_angle_surf_area (struct All_variables *E)
{
int es,el,m,i,j,ii,ia[5],lev;
- double aa,y1[4],y2[4],angle[6],xx[4][5],area_sphere_cap();
- void get_angle_sphere_cap();
+ double aa,y1[4],y2[4],angle[6],xx[4][5];
void parallel_process_termination();
for (m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -123,12 +119,10 @@
/* ================================================
area of spherical rectangle
================================================ */
-double area_sphere_cap(angle)
- double angle[6];
+double area_sphere_cap(double angle[6])
{
double area,a,b,c;
- double area_of_sphere_triag();
a = angle[1];
b = angle[2];
@@ -146,8 +140,7 @@
/* ================================================
area of spherical triangle
================================================ */
-double area_of_sphere_triag(a,b,c)
- double a,b,c;
+double area_of_sphere_triag(double a, double b, double c)
{
double ss,ak,aa,bb,cc,area;
@@ -176,13 +169,10 @@
angle [i]: angle bet test node and node i of the rectangle
angle1[i]: angle bet nodes i and i+1 of the rectangle
====================================================================== */
-double area_of_5points(E,lev,m,el,x,ne)
- struct All_variables *E;
- int lev,m,el,ne;
- double x[4];
+double area_of_5points(struct All_variables *E, int lev, int m, int el, double x[4], int ne)
{
int i,es,ia[5];
- double area1,get_angle(),area_of_sphere_triag();
+ double area1;
double xx[4],angle[5],angle1[5];
for (i=1;i<=4;i++)
@@ -210,12 +200,11 @@
get the angle for given four points spherical rectangle
================================= */
-void get_angle_sphere_cap(xx,angle)
- double xx[4][5],angle[6];
+void get_angle_sphere_cap(double xx[4][5], double angle[6])
{
int i,j,ii;
- double y1[4],y2[4],get_angle();;
+ double y1[4],y2[4];
for (i=1;i<=4;i++) { /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
for (j=1;j<=3;j++) {
@@ -238,8 +227,7 @@
/* ================================
get the angle for given two points
================================= */
-double get_angle(x,xx)
- double x[4],xx[4];
+double get_angle(double x[4], double xx[4])
{
double dist,angle;
const double pt5 = 0.5;
Modified: mc/3D/CitcomS/branches/eheien/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Stokes_flow_Incomp.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Stokes_flow_Incomp.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -57,9 +57,7 @@
/* Master loop for pressure and (hence) velocity field */
-void solve_constrained_flow_iterative(E)
- struct All_variables *E;
-
+void solve_constrained_flow_iterative(struct All_variables *E)
{
void v_from_vector();
void v_from_vector_pseudo_surf();
@@ -103,8 +101,8 @@
const int neq = E->lmesh.neq;
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- r1[m] = malloc((neq+1)*sizeof(double));
- r2[m] = malloc((neq+1)*sizeof(double));
+ r1[m] = (double*)malloc((neq+1)*sizeof(double));
+ r2[m] = (double*)malloc((neq+1)*sizeof(double));
}
/* r2 = F - grad(P) - K*V */
Modified: mc/3D/CitcomS/branches/eheien/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Topo_gravity.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Topo_gravity.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -63,11 +63,7 @@
*/
-void get_STD_topo(E,tpg,tpgb,divg,vort,ii)
- struct All_variables *E;
- float **tpg,**tpgb;
- float **divg,**vort;
- int ii;
+void get_STD_topo(struct All_variables *E, float **tpg, float **tpgb, float **divg, float **vort, int ii)
{
void allocate_STD_mem();
void compute_nodal_stress();
@@ -933,8 +929,7 @@
-void compute_geoid(E)
- struct All_variables *E;
+void compute_geoid(struct All_variables *E)
{
int i, p;
@@ -987,10 +982,7 @@
separately if stress output is needed
*/
-void get_CBF_topo(E,H,HB) /* call this only for top and bottom processors*/
- struct All_variables *E;
- float **H,**HB;
-
+void get_CBF_topo(struct All_variables *E, float **H, float **HB)
{
void get_elt_k();
void get_elt_g();
Modified: mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/Viscosity_structures.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -279,10 +279,7 @@
visc = E->VI[E->mesh.levmax]
*/
-void get_system_viscosity(E,propogate,evisc,visc)
- struct All_variables *E;
- int propogate;
- float **evisc,**visc;
+void get_system_viscosity(struct All_variables *E, int propogate, float **evisc, float **visc)
{
void visc_from_mat();
void visc_from_T();
@@ -408,9 +405,7 @@
}
-void visc_from_mat(E,EEta)
- struct All_variables *E;
- float **EEta;
+void visc_from_mat(struct All_variables *E, float **EEta)
{
int i,m,jj;
@@ -476,10 +471,7 @@
}
-void visc_from_T(E,EEta,propogate)
- struct All_variables *E;
- float **EEta;
- int propogate;
+void visc_from_T(struct All_variables *E, float **EEta, int propogate)
{
int m,i,k,l,z,jj,kk;
float zero,one,eta0,temp,tempa,TT[9];
@@ -575,8 +567,8 @@
computation of
depth, not needed
TWB */
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ TT[kk]=fmax(TT[kk],zero);
+ temp += fmin(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
}
EEta[m][ (i-1)*vpts + jj ] = tempa*
exp( E->viscosity.E[l]/(temp+E->viscosity.T[l])
@@ -604,8 +596,8 @@
temp=0.0;
zzz=0.0;
for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ TT[kk]=fmax(TT[kk],zero);
+ temp += fmin(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
}
@@ -635,8 +627,8 @@
for(jj=1;jj<=vpts;jj++) {
temp=0.0;
for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ TT[kk]=fmax(TT[kk],zero);
+ temp += fmin(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
}
if(E->control.mat_control==0){
@@ -686,8 +678,8 @@
for(jj=1;jj <= vpts;jj++) {
temp=0.0;zzz=0.0;
for(kk=1;kk <= ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ TT[kk]=fmax(TT[kk],zero);
+ temp += fmin(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
}
EEta[m][ (i-1)*vpts + jj ] = tempa*
@@ -797,8 +789,8 @@
for(jj=1;jj<=vpts;jj++) {
temp=zzz=0.0;
for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)]; /* mean temp */
+ TT[kk]=fmax(TT[kk],zero);
+ temp += fmin(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)]; /* mean temp */
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];/* mean r */
}
/* convert to z, as defined to be unity at surface
@@ -902,10 +894,7 @@
}
-void visc_from_S(E,EEta,propogate)
- struct All_variables *E;
- float **EEta;
- int propogate;
+void visc_from_S(struct All_variables *E, float **EEta, int propogate)
{
float one,two,scale,stress_magnitude,depth,exponent1;
float *eedot;
@@ -933,7 +922,7 @@
}
/* eedot cannot be too small, or the viscosity will go to inf */
for(e=1;e<=nel;e++){
- eedot[e] = max(eedot[e], 1.0e-16);
+ eedot[e] = fmax(eedot[e], 1.0e-16);
}
for(e=1;e<=nel;e++) {
@@ -948,7 +937,7 @@
return;
}
-void visc_from_P(E,EEta) /* "plasticity" implementation
+void visc_from_P(struct All_variables *E, float **EEta) /* "plasticity" implementation
psrw = FALSE
@@ -993,8 +982,6 @@
TWB
*/
- struct All_variables *E;
- float **EEta;
{
float *eedot,zz[9],zzz,tau,eta_p,eta_new,tau2,eta_old,eta_old2;
int m,e,l,z,jj,kk;
@@ -1047,7 +1034,7 @@
tau = E->viscosity.pdepv_a[l] + zzz * E->viscosity.pdepv_b[l];
/* min of depth dep. and constant yield stress */
- tau = min(tau, E->viscosity.pdepv_y[l]);
+ tau = fmin(tau, E->viscosity.pdepv_y[l]);
/* yield viscosity */
eta_p = tau/(2.0 * eedot[e] + 1e-7) + E->viscosity.pdepv_offset;
@@ -1056,7 +1043,7 @@
eta_new = 1.0/(1.0/EEta[m][ (e-1)*vpts + jj ] + 1.0/eta_p);
}else{
/* min viscosities*/
- eta_new = min(EEta[m][ (e-1)*vpts + jj ], eta_p);
+ eta_new = fmin(EEta[m][ (e-1)*vpts + jj ], eta_p);
}
//fprintf(stderr,"z: %11g mat: %i a: %11g b: %11g y: %11g ee: %11g tau: %11g eta_p: %11g eta_new: %11g eta_old: %11g\n",
// zzz,l,E->viscosity.pdepv_a[l], E->viscosity.pdepv_b[l],E->viscosity.pdepv_y[l],
@@ -1077,7 +1064,7 @@
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
/* compute sigma_y as above */
tau = E->viscosity.pdepv_a[l] + zzz * E->viscosity.pdepv_b[l];
- tau = min(tau, E->viscosity.pdepv_y[l]);
+ tau = fmin(tau, E->viscosity.pdepv_y[l]);
tau2 = tau * tau;
if(tau < 1e10){
/* */
@@ -1105,9 +1092,7 @@
compositions between zero and unity
*/
-void visc_from_C( E, EEta)
- struct All_variables *E;
- float **EEta;
+void visc_from_C(struct All_variables *E, float **EEta)
{
double vmean,cc_loc[10],CC[10][9],cbackground;
int m,l,z,jj,kk,i,p,q;
@@ -1155,10 +1140,7 @@
} /* end cap */
}
-void strain_rate_2_inv(E,m,EEDOT,SQRT)
- struct All_variables *E;
- float *EEDOT;
- int m,SQRT;
+void strain_rate_2_inv(struct All_variables *E, int m, float *EEDOT, int SQRT)
{
void get_rtf_at_ppts();
void velo_from_element();
@@ -1397,7 +1379,7 @@
E->sx[m][3][E->ien[m][ee].node[8]]);
/* if ee has tracers in it and is within the channel */
- if((E->trace.ntracer_flavor[m][flavor][ee] > 0) &&
+ if((E->trace.num_tracer_flavors[m][flavor][ee] > 0) &&
(rad_mean <= rr + E->viscosity.lv_channel_thickness)) {
F[e] = E->viscosity.lv_reduction;
break;
@@ -1458,7 +1440,7 @@
ee = (k-1)*E->lmesh.elz + ii;
/* if ee has tracers in it */
- if(E->trace.ntracer_flavor[m][flavor][ee] > 0) {
+ if(E->trace.num_tracer_flavors[m][flavor][ee] > 0) {
F[e] = E->viscosity.lv_reduction;
break;
}
Modified: mc/3D/CitcomS/branches/eheien/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/global_defs.h 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/global_defs.h 2011-07-14 00:06:30 UTC (rev 18754)
@@ -51,23 +51,17 @@
#include "hdf5.h"
#endif
-#ifdef __cplusplus
+//#ifdef __cplusplus
+//extern "C" {
+//#endif
-extern "C" {
-#else
-
-
-
-
/* Macros */
-#define max(A,B) (((A) > (B)) ? (A) : (B))
-#define min(A,B) (((A) < (B)) ? (A) : (B))
+//#define max(A,B) (((A) > (B)) ? (A) : (B))
+//#define min(A,B) (((A) < (B)) ? (A) : (B))
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
-#endif
-
#define LIDN 0x1
#define VBX 0x2
#define VBZ 0x4
@@ -700,8 +694,17 @@
#ifdef USE_HDF5
#include "hdf5_related.h"
#endif
+
+//#ifdef __cplusplus
+//}
+//#endif
+
#include "tracer_defs.h"
+//#ifdef __cplusplus
+//extern "C" {
+//#endif
+
struct All_variables {
#include "solver.h"
@@ -855,9 +858,9 @@
};
-#ifdef __cplusplus
-}
-#endif
+//#ifdef __cplusplus
+//}
+//#endif
#ifndef TRUE
#define TRUE 1
Modified: mc/3D/CitcomS/branches/eheien/lib/prototypes.h
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/prototypes.h 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/prototypes.h 2011-07-14 00:06:30 UTC (rev 18754)
@@ -1,3 +1,7 @@
+#ifdef __cplusplus
+extern "C" {
+#endif
+
/* Advection_diffusion.c */
void advection_diffusion_parameters(struct All_variables *);
void advection_diffusion_allocate_memory(struct All_variables *);
@@ -528,3 +532,8 @@
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]);
+
+#ifdef __cplusplus
+}
+#endif
+
Modified: mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/lib/tracer_defs.h 2011-07-14 00:06:30 UTC (rev 18754)
@@ -26,12 +26,32 @@
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
+#include <vector>
+#include <list>
-
/* forward declaration */
struct All_variables;
+class Tracer {
+public:
+ // Tracer position in spherical coordinates
+ double theta, phi, rad;
+ // Tracer position in Cartesian coordinates
+ double x, y, z;
+ // Previous Cartesian position and velocity
+ double x0, y0, z0;
+ double Vx, Vy, Vz;
+
+ // Tracer flavor (meaning should be application dependent)
+ double flavor;
+
+ int ielement;
+};
+
+typedef std::list<Tracer> TracerList;
+typedef std::vector<TracerList> TracerArray;
+
struct TRACE{
FILE *fpt;
@@ -47,26 +67,27 @@
double box_cushion;
/* tracer arrays */
- int number_of_basic_quantities;
- int number_of_extra_quantities;
- int number_of_tracer_quantities;
+ TracerArray tracers;
+ //int number_of_basic_quantities;
+ //int number_of_extra_quantities;
+ //int number_of_tracer_quantities;
- double *basicq[13][100];
- double *extraq[13][100];
+ //double *basicq[13][100];
+ //double *extraq[13][100];
- int ntracers[13];
- int max_ntracers[13];
- int *ielement[13];
+ //int ntracers[13];
+ //int max_ntracers[13];
+ //int *ielement[13];
- int number_of_tracers;
+ //int ilatersize[13];
+ //int ilater[13];
+ //double *rlater[13][100];
+ TracerList later_tracer[13];
- int ilatersize[13];
- int ilater[13];
- double *rlater[13][100];
-
/* tracer flavors */
int nflavors;
- int **ntracer_flavor[13];
+ std::map<int, std::map<int, int> > num_tracer_flavors[13];
+ //int **ntracer_flavor[13];
int ic_method_for_flavors;
double *z_interface;
Modified: mc/3D/CitcomS/branches/eheien/module/misc.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/module/misc.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/module/misc.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -123,7 +123,7 @@
"%s: 'libCitcomSCommon.citcom_init' failed",
pyCitcom_citcom_init__name__);
- cobj = PyCObject_FromVoidPtr(E, deleteE);
+ cobj = PyCObject_FromVoidPtr(E, (void (*)(void*))deleteE);
return Py_BuildValue("N", cobj);
}
Modified: mc/3D/CitcomS/branches/eheien/visual/project_geoid.c
===================================================================
--- mc/3D/CitcomS/branches/eheien/visual/project_geoid.c 2011-07-13 19:52:44 UTC (rev 18753)
+++ mc/3D/CitcomS/branches/eheien/visual/project_geoid.c 2011-07-14 00:06:30 UTC (rev 18754)
@@ -66,10 +66,10 @@
void allocate_sph_harm(sph_harm *coeff, int len)
{
- coeff->ll = malloc(len * sizeof(int));
- coeff->mm = malloc(len * sizeof(int));
- coeff->clm = malloc(len * sizeof(float));
- coeff->slm = malloc(len * sizeof(float));
+ coeff->ll = (int*)malloc(len * sizeof(int));
+ coeff->mm = (int*)malloc(len * sizeof(int));
+ coeff->clm = (float*)malloc(len * sizeof(float));
+ coeff->slm = (float*)malloc(len * sizeof(float));
}
@@ -145,8 +145,8 @@
grid->nphi = nphi;
/* allocate memory */
- grid->theta = malloc(ntheta * sizeof(double));
- grid->phi = malloc(nphi * sizeof(double));
+ grid->theta = (double*)malloc(ntheta * sizeof(double));
+ grid->phi = (double*)malloc(nphi * sizeof(double));
/* create a regular mesh */
More information about the CIG-COMMITS
mailing list