[cig-commits] r18261 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Thu Apr 21 06:04:39 PDT 2011
Author: becker
Date: 2011-04-21 06:04:38 -0700 (Thu, 21 Apr 2011)
New Revision: 18261
Modified:
mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
mc/3D/CitcomCU/trunk/src/Citcom.c
mc/3D/CitcomCU/trunk/src/Construct_arrays.c
mc/3D/CitcomCU/trunk/src/Drive_solvers.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Makefile
mc/3D/CitcomCU/trunk/src/Output_gzdir.c
mc/3D/CitcomCU/trunk/src/Parallel_related.c
mc/3D/CitcomCU/trunk/src/Solver_multigrid.c
mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
mc/3D/CitcomCU/trunk/src/global_defs.h
mc/3D/CitcomCU/trunk/src/prototypes.h
mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
Log:
Experimental change of surface BC from free slip to no slip (not working
yet, but checked in for now), and mode to assign compositional viscosity
as an absolute value rather than factor
Modified: mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -121,6 +121,30 @@
return;
}
+
+
+void freeze_surface(struct All_variables *E)
+{
+
+ int lev,top;
+ if(E->parallel.me == 0)
+ fprintf(stderr,"WARNING: freezing surface boundary condition at time step %i\n",
+ E->monitor.solution_cycles);
+ /* no slip on top */
+ E->mesh.topvbc = 1;E->control.VBXtopval=E->control.VBYtopval=E->control.plate_vel=0.0;
+ velocity_boundary_conditions(E);
+
+ check_bc_consistency(E);
+ construct_id(E);
+ //construct_lm(E);
+ //construct_sub_element(E);
+ parallel_shuffle_ele_and_id(E);
+ parallel_communication_routs(E);
+ //construct_shape_functions(E);
+ //mass_matrix(E);
+ velocities_conform_bcs(E, E->U);
+
+}
/* ========================================== */
void temperature_boundary_conditions(struct All_variables *E)
Modified: mc/3D/CitcomCU/trunk/src/Citcom.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Citcom.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Citcom.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -110,30 +110,38 @@
while(E.control.keep_going && (Emergency_stop == 0))
{
- //if(E.parallel.me == 0)fprintf(stderr,"processing heating\n");
+
+
process_heating(&E);
+ if(E.parallel.me==0)fprintf(stderr,"process heating done\n");
+
E.monitor.solution_cycles++;
if(E.monitor.solution_cycles > E.control.print_convergence)
E.control.print_convergence = 1;
- //if(E.parallel.me == 0)fprintf(stderr,"processing buoyancy\n");
/**/ report(&E, "Update buoyancy for further `timesteps'");
(E.next_buoyancy_field) (&E);
- //if(E.parallel.me == 0)fprintf(stderr,"processing temp\n");
+ if(E.parallel.me==0)fprintf(stderr,"buoyancy field done\n");
+
/**/ report(&E, "Process results of buoyancy update");
process_temp_field(&E, E.monitor.solution_cycles);
+ if(E.parallel.me==0)fprintf(stderr,"temp field done\n");
- //if(E.parallel.me == 0)fprintf(stderr,"solving stokes\n");
-
+ if(E.monitor.solution_cycles == E.control.freeze_surface_at_step){ /* for testing purposes */
+ freeze_surface(&E);
+ }
+
general_stokes_solver(&E);
-
+ if(E.parallel.me==0)fprintf(stderr,"stokes solver done\n");
+
if(E.control.composition){
- //if(E.parallel.me == 0)fprintf(stderr,"composition, update buoyancy with R-G\n");
(E.next_buoyancy_field) (&E); /* correct with R-G */
+ if(E.parallel.me==0)fprintf(stderr,"next buoyancy composition done\n");
}
- //if(E.parallel.me == 0)fprintf(stderr,"processing new velocity\n");
/**/ report(&E, "Process results of velocity solver");
process_new_velocity(&E, E.monitor.solution_cycles);
+ if(E.parallel.me==0)fprintf(stderr,"process new velocity done\n");
+
if(E.monitor.T_interior > E.monitor.T_interior_max)
{
fprintf(stderr, "quit due to maxT = %.4e sub_iteration%d\n", E.monitor.T_interior, E.advection.last_sub_iterations);
Modified: mc/3D/CitcomCU/trunk/src/Construct_arrays.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Construct_arrays.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -144,8 +144,6 @@
int elx, ely, elz;
const int dims = E->mesh.nsd;
- //const int dofs = E->mesh.dof;
- //const int ends = enodes[dims];
for(lev = E->mesh.levmax; lev >= E->mesh.levmin; lev--)
{
@@ -165,7 +163,7 @@
eqn_count += 1;
}
}
-
+
E->lmesh.NEQ[lev] = eqn_count;
j1 = 1;
@@ -192,11 +190,11 @@
E->parallel.IDD[lev][E->ID[lev][node].doff[doff]] = 1;
}
+
}
E->lmesh.neq = E->lmesh.NEQ[E->mesh.levmax]; /* Total NUMBER of independent variables */
-
lev = E->mesh.levmax;
if(E->control.verbose)
{
@@ -282,7 +280,7 @@
void construct_node_maps(struct All_variables *E)
{
//float initial_time;
-
+ static int been_here = 0;
//int el, n, nn, lev, i, j, k, ja, jj, ii, kk, ia, is, ie, js, je, ks, ke, dims2;
int nn, lev, i, j, k, ja, jj, ii, kk, ia, is, ie, js, je, ks, ke;
//int doff, nox, noy, noz, noxz, node1, eqn1, loc1, count, found, element;
@@ -342,15 +340,18 @@
}
}
}
+ if(!been_here){ /* allow having this routine called more than once */
+ E->Eqn_k1[lev] = (higher_precision *) malloc((matrix + 5) * sizeof(higher_precision));
+ E->Eqn_k2[lev] = (higher_precision *) malloc((matrix + 5) * sizeof(higher_precision));
+ if(dims == 3)
+ E->Eqn_k3[lev] = (higher_precision *) malloc((matrix + 5) * sizeof(higher_precision));
- E->Eqn_k1[lev] = (higher_precision *) malloc((matrix + 5) * sizeof(higher_precision));
- E->Eqn_k2[lev] = (higher_precision *) malloc((matrix + 5) * sizeof(higher_precision));
- if(dims == 3)
- E->Eqn_k3[lev] = (higher_precision *) malloc((matrix + 5) * sizeof(higher_precision));
+ }
E->mesh.matrix_size[lev] = matrix + 1;
} /* end for level and m */
+ been_here = 1;
return;
}
@@ -364,7 +365,6 @@
int neq, nno, nel;
double elt_K[24 * 24];
- //static int been_here = 0;
double w1, w2, w3, ww1, ww2, ww3;
//higher_precision *B1, *B2, *B3;
@@ -846,7 +846,8 @@
E->elt_k[i] = (struct EK *)malloc((E->lmesh.NEL[i] + 1) * sizeof(struct EK));
}
}
- if(been_here0 == 0 || (E->viscosity.update_allowed && E->monitor.solution_cycles % E->control.KERNEL == 0))
+ if((been_here0 == 0) || (E->viscosity.update_allowed && E->monitor.solution_cycles % E->control.KERNEL == 0) ||
+ (E->monitor.solution_cycles == E->control.freeze_surface_at_step) )
{
/* do the following for the 1st time or update_allowed is true */
@@ -857,8 +858,10 @@
if(E->control.NMULTIGRID || E->control.NASSEMBLE)
{
- if(been_here == 0)
+ if((been_here == 0)||(E->monitor.solution_cycles == E->control.freeze_surface_at_step) )
{ /* node_maps only built once */
+ if((E->parallel.me==0) && (E->monitor.solution_cycles == E->control.freeze_surface_at_step))
+ fprintf(stderr,"reconstructing node maps\n");
construct_node_maps(E);
been_here = 1;
}
Modified: mc/3D/CitcomCU/trunk/src/Drive_solvers.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Drive_solvers.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -111,7 +111,7 @@
get_system_viscosity(E, 1, E->EVI[E->mesh.levmax], E->VI[E->mesh.levmax]);
construct_stiffness_B_matrix(E);
- //if(E->parallel.me==0)fprintf(stderr,"calling solver\n");
+ if(E->parallel.me==0)fprintf(stderr,"calling solver\n");
solve_constrained_flow_iterative(E);
@@ -146,7 +146,7 @@
if(iterate){ /* free the delta_U array */
free((void *)delta_U);
}
- //if(E->parallel.me==0)fprintf(stderr,"stokes solver done\n");
+ if(E->parallel.me==0)fprintf(stderr,"stokes solver done\n");
return;
}
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -165,6 +165,7 @@
}
+
/* ===================================
Functions which set up details
common to all problems follow ...
@@ -809,6 +810,8 @@
input_int("stokes_flow_only", &(E->control.stokes), "0", m);
+ input_int("freeze_surface_at_step",&(E->control.freeze_surface_at_step),"-9999999",m);
+
input_boolean("node_assemble", &(E->control.NASSEMBLE), "off", m);
/* general mesh structure */
Modified: mc/3D/CitcomCU/trunk/src/Makefile
===================================================================
--- mc/3D/CitcomCU/trunk/src/Makefile 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Makefile 2011-04-21 13:04:38 UTC (rev 18261)
@@ -107,7 +107,8 @@
LinuxFLAGS=
LinuxLDFLAGS=
#LinuxOPTIM=-g
-LinuxOPTIM=-O2
+#LinuxOPTIM=-O2
+LinuxOPTIM=-O2 -Wall -Wextra -pedantic -std=c99 -Wunused-function -D_GNU_SOURCE
####################################
#PARAGON
@@ -145,8 +146,11 @@
Global_operations.c\
Stokes_flow_Incomp.c\
Instructions.c\
+ io.c \
Nodal_mesh.c\
Output.c\
+ output_ascii.c \
+ output_vtk.c \
Pan_problem_misc_functions.c\
Parallel_related.c\
Parsing.c\
@@ -241,9 +245,18 @@
Instructions.o: $(HEADER) Instructions.c
$(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Instructions.c
+io.o: $(HEADER) io.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) io.c
+
Nodal_mesh.o: $(HEADER) Nodal_mesh.c
$(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Nodal_mesh.c
+output_ascii.o: $(HEADER) output_ascii.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) output_ascii.c
+
+output_vtk.o: $(HEADER) output_vtk.c
+ $(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) output_vtk.c
+
Output.o: $(HEADER) Output.c
$(CC) $(OPTIM) $(FLAGS) $(OBJFLAG) Output.c
Modified: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -651,7 +651,7 @@
}
}
- //if(E->parallel.me==0)fprintf(stderr,"vel output done\n");
+ if(E->parallel.me==0)fprintf(stderr,"vel output done\n");
parallel_process_sync();
return;
Modified: mc/3D/CitcomCU/trunk/src/Parallel_related.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Parallel_related.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Parallel_related.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -1009,9 +1009,6 @@
//int target_proc, kk, e, node, i, ii, j, k, bound, type, idb, msginfo[8];
int target_proc, kk, k, idb;
- //static int been_here = 0;
- //static int mid_recv, sizeofk;
- //const int levmax = E->mesh.levmax;
MPI_Status status[100];
MPI_Request request[100];
Modified: mc/3D/CitcomCU/trunk/src/Solver_multigrid.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Solver_multigrid.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Solver_multigrid.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -312,7 +312,6 @@
//float average, w;
float average;
- //static int been_here = 0;
const int sl_minus = start_lev - 1;
const int nels_minus = E->lmesh.NEL[start_lev - 1];
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2011-04-21 13:04:38 UTC (rev 18261)
@@ -60,7 +60,7 @@
int m = E->parallel.me ;
/* default values .... */
E->viscosity.update_allowed = 0;
- E->viscosity.SDEPV = E->viscosity.TDEPV = E->viscosity.BDEPV =
+ E->viscosity.SDEPV = E->viscosity.TDEPV = E->viscosity.BDEPV =
E->viscosity.CDEPV = 0;
E->viscosity.EXPX = 0;
E->viscosity.SMOOTH = 0;
@@ -149,8 +149,14 @@
input_boolean("TDEPV", &(E->viscosity.TDEPV), "on",m);
input_boolean("SDEPV", &(E->viscosity.SDEPV), "off",m);
+ input_int("sdepv_rheology", &(E->viscosity.sdepv_rheology), "1",m); /* type of stress dependence */
+
+ input_boolean("sdepv_start_from_newtonian", &(E->viscosity.sdepv_start_from_newtonian), "off",m);
+
+
input_boolean("BDEPV",&(E->viscosity.BDEPV),"off",m);
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
+ input_boolean("cdepv_absolute",&(E->viscosity.cdepv_absolute),"off",m); /* make comp dep viscosity absolute */
/* plasticity offset viscosity */
input_float("plasticity_viscosity_offset", &(E->viscosity.plasticity_viscosity_offset),"0.0",m);
@@ -330,6 +336,7 @@
if(E->viscosity.CDEPV)
visc_from_C(E, visc, evisc, propogate);
+ /* two different kinds of implementations of powerlaw */
if(E->viscosity.SDEPV)
visc_from_S(E, visc, evisc, propogate);
@@ -693,16 +700,34 @@
}
else
strain_rate_2_inv(E, eedot, 1);
-
- for(e = 1; e <= nel; e++)
- {
+ if((!E->viscosity.sdepv_start_from_newtonian)||(visits)){
+ switch(E->viscosity.sdepv_rheology){
+ case 1: /* old default, i think the factors
+ don't make sense, leave in for
+ backward compatibility */
+ for(e = 1; e <= nel; e++)
+ {
exponent1 = one - one / E->viscosity.sdepv_expt[E->mat[e] - 1];
scale = pow(two * eedot[e] / E->viscosity.sdepv_trns[E->mat[e] - 1], exponent1);
for(jj = 1; jj <= vpts; jj++)
- EEta[(e - 1) * vpts + jj] = two * EEta[(e - 1) * vpts + jj] / (one + scale * pow(EEta[(e - 1) * vpts + jj], exponent1));
+ EEta[(e - 1) * vpts + jj] = two * EEta[(e - 1) * vpts + jj] / (one + scale * pow(EEta[(e - 1) * vpts + jj], exponent1));
+ }
+ break;
+ case 2: /* composite rheology, different from
+ above by the missing two up front*/
+ for(e = 1; e <= nel; e++){
+ exponent1 = one - one / E->viscosity.sdepv_expt[E->mat[e] - 1];
+ scale = pow(two * eedot[e] / E->viscosity.sdepv_trns[E->mat[e] - 1], exponent1);
+ for(jj = 1; jj <= vpts; jj++)
+ EEta[(e - 1) * vpts + jj] = EEta[(e - 1) * vpts + jj] / (one + scale * pow(EEta[(e - 1) * vpts + jj], exponent1));
+ }
+
+ break;
+ default:
+ myerror("stress dependent rheology mode undefined",E);
+ }
}
-
visits++;
free((void *)eedot);
@@ -1409,8 +1434,37 @@
/* log of the material viscosities */
logv[0] = log(E->viscosity.pre_comp[0]);
logv[1] = log(E->viscosity.pre_comp[1]);
+ if((E->parallel.me==0)&&(E->viscosity.cdepv_absolute)){
+ fprintf(stderr,"WARNING: using cdepv for absolute viscosity\n");
+ }
}
- for(i = 1; i <= nel; i++)
+ if(E->viscosity.cdepv_absolute){
+ /*
+
+ override all other viscosities (this is exact same as below but
+ for the viscosity assignment (repeated here to avoid if statments)
+ */
+ for(i = 1; i <= nel; i++){
+ for(kk = 1; kk <= ends; kk++){
+ CC[kk] = E->C[E->ien[i].node[kk]];
+ if(E->control.check_c_irange){
+ if(CC[kk] < 0)CC[kk]=0.0;if(CC[kk] > 1)CC[kk]=1.0;
+ }
+ }
+ for(jj = 1; jj <= vpts; jj++){
+ cc_loc = 0.0;
+ for(kk = 1; kk <= ends; kk++)
+ cc_loc += CC[kk] * E->N.vpt[GNVINDEX(kk, jj)];
+ vmean = exp(cc_loc * logv[1] + (1.0-cc_loc) * logv[0]);
+ EEta[ (i-1)*vpts + jj ] = vmean;
+ } /* end jj loop */
+ }
+ }else{
+ /*
+ regular mode, multiply with any previous viscosity
+
+ */
+ for(i = 1; i <= nel; i++)
{
/* determine composition of each of the nodes of the element */
for(kk = 1; kk <= ends; kk++){
@@ -1432,6 +1486,7 @@
EEta[ (i-1)*vpts + jj ] *= vmean;
} /* end jj loop */
}
+ }
visited++;
}
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2011-04-21 13:04:38 UTC (rev 18261)
@@ -750,8 +750,8 @@
int CONMAN;
int stokes;
+ int freeze_surface_at_step;
-
int check_t_irange,check_c_irange;
float Ra_670, clapeyron670, transT670, width670;
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2011-04-21 13:04:38 UTC (rev 18261)
@@ -40,6 +40,7 @@
void drex_eigen(double [3][3], double [3][3], int *);
void malmul_scaled_id(double [3][3], double [3][3], double, double);
/* Boundary_conditions.c */
+void freeze_surface(struct All_variables *);
void velocity_boundary_conditions(struct All_variables *);
void temperature_boundary_conditions(struct All_variables *);
void velocity_refl_vert_bc(struct All_variables *);
@@ -361,6 +362,7 @@
void visc_from_mat(struct All_variables *, float *, float *);
void visc_from_T(struct All_variables *, float *, float *, int);
void visc_from_S(struct All_variables *, float *, float *, int);
+void visc_from_S2(struct All_variables *, float *, float *, int);
void strain_rate_2_inv(struct All_variables *, float *, int);
double second_invariant_from_3x3(double [3][3]);
void calc_vgm_matrix(struct All_variables *, double *, double *);
Modified: mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2011-04-20 22:51:03 UTC (rev 18260)
+++ mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2011-04-21 13:04:38 UTC (rev 18261)
@@ -93,10 +93,12 @@
int SDEPV;
float sdepv_misfit;
+ int sdepv_rheology;
float sdepv_iter_damp;
int sdepv_normalize;
float sdepv_expt[CITCOM_CU_VISC_MAXLAYER];
float sdepv_trns[CITCOM_CU_VISC_MAXLAYER];
+ int sdepv_start_from_newtonian;
int TDEPV;
int TDEPV_AVE;
@@ -139,6 +141,8 @@
int CDEPV; /* composition dependent viscosity */
+ int cdepv_absolute; /* make the composition an absolute viscosity,
+ not a prefactor */
float pre_comp[CITCOM_CU_VISC_MAXLAYER]; /* prefactors */
More information about the CIG-COMMITS
mailing list