[cig-commits] r19936 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Tue Apr 10 16:55:53 PDT 2012
Author: becker
Date: 2012-04-10 16:55:52 -0700 (Tue, 10 Apr 2012)
New Revision: 19936
Modified:
mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Instructions.c
mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
mc/3D/CitcomCU/trunk/src/Process_velocity.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:
Composition dependent plasticity
Testing implementation of periodic boundary conditions now includes a net drift
filter (still testing)
Modified: mc/3D/CitcomCU/trunk/src/Boundary_conditions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/Boundary_conditions.c 2012-04-10 23:55:52 UTC (rev 19936)
@@ -70,9 +70,11 @@
}
}
- if(E->mesh.periodic_x || E->mesh.periodic_y)
+ if(E->mesh.periodic_x || E->mesh.periodic_y){
+
velocity_apply_periodic_bcs(E);
- else
+
+ }else
velocity_refl_vert_bc(E); /* default */
for(lv = E->mesh.levmax; lv >= E->mesh.levmin; lv--)
@@ -99,7 +101,31 @@
}
}
-
+ if(E->mesh.periodic_pin_or_filter == 1){
+ if(E->mesh.periodic_x || E->mesh.periodic_y){
+
+
+ /*
+ pin one node at top in lower left corner
+ */
+ for(lv = E->mesh.levmax; lv >= E->mesh.levmin; lv--){
+ node = E->mesh.NOZ[lv];
+ if(E->mesh.periodic_x){
+ E->NODE[lv][node] = E->NODE[lv][node] & (~SBX); /* turn of stress */
+ E->NODE[lv][node] = E->NODE[lv][node] | (VBX);/* turn on velocity BC */
+ if(lv == E->mesh.levmax){ /* NB */
+ E->VB[1][node] = 0; /* set to zero */
+ }
+ }else{ /* assume only x or y periodic */
+ E->NODE[lv][node] = E->NODE[lv][node] & (~SBY);
+ E->NODE[lv][node] = E->NODE[lv][node] | (VBY);
+ if(lv == E->mesh.levmax){ /* NB */
+ E->VB[2][node] = 0;
+ }
+ }
+ }
+ }
+ }
if(E->control.verbose)
{
for(node = 1; node <= E->lmesh.nno; node++)
@@ -467,11 +493,11 @@
fprintf(E->fp,"Periodic boundary conditions\n");
- if (E->mesh.periodic_y && E->mesh.periodic_x) {
- return;
- }
+ //if (E->mesh.periodic_y && E->mesh.periodic_x) {
+ // return;
+ //}
- else if (E->mesh.periodic_y) {
+ if (E->mesh.periodic_y) {
/* except one side with XOZ and y=0, all others are not reflecting BC*/
/* for two YOZ planes if 3-D, or two OZ side walls for 2-D */
@@ -497,7 +523,7 @@
}
- else if (E->mesh.periodic_x) {
+ if (E->mesh.periodic_x) {
/* for two XOZ planes if 3-D */
@@ -611,11 +637,11 @@
/* Temps and bc-values at top level only */
/* fixed temperature at x=0 */
- if (E->mesh.periodic_x && E->mesh.periodic_y) {
- return;
- }
+ //if (E->mesh.periodic_x && E->mesh.periodic_y) {
+ //return;
+ //}
- else if (E->mesh.periodic_y) {
+ if (E->mesh.periodic_y) {
if (E->parallel.me_loc[1]==0 || E->parallel.me_loc[1]==E->parallel.nprocx-1)
for(j=1;j<=E->lmesh.noy;j++)
@@ -635,7 +661,7 @@
} /* end for loop i & j */
}
- else if (E->mesh.periodic_x) {
+ if (E->mesh.periodic_x) {
if (E->parallel.me_loc[2]==0 || E->parallel.me_loc[2]==E->parallel.nprocy-1)
for(j=1;j<=E->lmesh.nox;j++)
for(i=1;i<=E->lmesh.noz;i++) {
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2012-04-10 23:55:52 UTC (rev 19936)
@@ -52,6 +52,10 @@
{
input_boolean("composition", &(E->control.composition), "0", E->parallel.me);
+
+ input_boolean("composition_neutralize_buoyancy", &(E->control.composition_neutralize_buoyancy), "0", E->parallel.me);
+ if(E->control.composition_neutralize_buoyancy && E->parallel.me == 0)
+ fprintf(stderr,"WARNING: reducing thermal buoyanyc to zero in composition=1 regions\n");
input_int("tracers_add_flavors", &(E->tracers_add_flavors), "0", E->parallel.me);
input_float("tracers_assign_dense_fraction",&(E->tracers_assign_dense_fraction),"1.0",E->parallel.me);
Modified: mc/3D/CitcomCU/trunk/src/Instructions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Instructions.c 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/Instructions.c 2012-04-10 23:55:52 UTC (rev 19936)
@@ -454,6 +454,7 @@
E->control.CHEMISTRY_MODULE = 0;
E->control.composition = 0;
+ E->control.composition_neutralize_buoyancy = 0;
E->sphere.vtk_base_init = 0;
@@ -495,6 +496,8 @@
E->mesh.sidevbc = 0;
E->mesh.periodic_x = 0; /* reflection is default */
E->mesh.periodic_y = 0;
+ E->mesh.periodic_pin_or_filter = 0; /* filter by default */
+
E->control.VBXtopval = 0.0;
E->control.VBYtopval = 0.0;
E->control.VBXbotval = 0.0;
@@ -927,6 +930,8 @@
input_boolean("periodicx", &(E->mesh.periodic_x), "off", m);
input_boolean("periodicy", &(E->mesh.periodic_y), "off", m);
+ input_int("periodic_pin_or_filter", &(E->mesh.periodic_pin_or_filter), "0", m); /* 1: pin 0: filter */
+
input_boolean("depthdominated", &(E->control.depth_dominated), "off", m);
input_boolean("eqnzigzag", &(E->control.eqn_zigzag), "off", m);
input_boolean("eqnviscosity", &(E->control.eqn_viscosity), "off", m);
Modified: mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/Pan_problem_misc_functions.c 2012-04-10 23:55:52 UTC (rev 19936)
@@ -85,10 +85,25 @@
//const int i2 = 670;
H = (float *)malloc((E->lmesh.noz + 1) * sizeof(float));
- for(i = 1; i <= E->lmesh.nno; i++)
- {
- j = (i - 1) % (E->lmesh.noz) + 1;
- E->buoyancy[i] = E->control.Atemp * E->T[i] * E->expansivity[j] - E->control.Acomp * E->C[i];
+ if(E->control.composition_neutralize_buoyancy){
+
+ for(i = 1; i <= E->lmesh.nno; i++){ /* for testing purposes,
+ have unity
+ composition reset the
+ buoyancy */
+ if(E->C[i]>1)E->C[i] = 1;
+ if(E->C[i]<0)E->C[i] = 0;
+ j = (i - 1) % (E->lmesh.noz) + 1;
+ E->buoyancy[i] = (1.-E->C[i]) * E->control.Atemp * E->T[i] * E->expansivity[j] ;
+
+ }
+
+ }else{ /* default */
+ for(i = 1; i <= E->lmesh.nno; i++)
+ {
+ j = (i - 1) % (E->lmesh.noz) + 1;
+ E->buoyancy[i] = E->control.Atemp * E->T[i] * E->expansivity[j] - E->control.Acomp * E->C[i];
+ }
}
if(E->control.Ra_670 != 0.0 || E->control.Ra_410 != 0.0)
{
Modified: mc/3D/CitcomCU/trunk/src/Process_velocity.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Process_velocity.c 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/Process_velocity.c 2012-04-10 23:55:52 UTC (rev 19936)
@@ -51,6 +51,14 @@
void process_new_velocity(struct All_variables *E, int ii)
{
+ if(E->mesh.periodic_pin_or_filter == 0){
+ /* check for rigid body motion */
+ if(E->mesh.periodic_x || E->mesh.periodic_y){
+ remove_net_drift(E);
+
+ }
+ }
+
if(E->control.stokes || ((ii % E->control.record_every) == 0))
{
/* get_CBF_topo(E,E->slice.tpg,E->slice.tpgb); */
@@ -76,6 +84,7 @@
}
}
+
return;
}
@@ -209,3 +218,28 @@
return;
}
+
+void remove_net_drift(struct All_variables *E)
+{
+ float net_drift;
+ int coord,i;
+
+ if(E->mesh.periodic_x && E->mesh.periodic_y)
+ myerror("remove_net_drift: can't deal with both x and y periodicity",E);
+
+ if(E->mesh.periodic_x)
+ coord = 1;
+ else if(E->mesh.periodic_y)
+ coord = 2;
+
+ net_drift = return_bulk_value(E, E->V[coord], -1,1);
+
+ if(E->parallel.me == 0)
+ fprintf(stderr,"remove_net_drift: removing net drift of %g in %i coord\n",
+ net_drift,coord);
+
+ for(i = 1; i <= E->lmesh.nno; i++)
+ E->V[coord][i] -= net_drift;
+
+}
+
Modified: mc/3D/CitcomCU/trunk/src/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/Viscosity_structures.c 2012-04-10 23:55:52 UTC (rev 19936)
@@ -49,6 +49,7 @@
#include "anisotropic_viscosity.h"
#endif
+
static void visc_from_B(struct All_variables *, float *, float *, int );
static void visc_from_C(struct All_variables *, float *, float *, int );
void *safe_malloc (size_t );
@@ -156,8 +157,11 @@
input_int("max_sdep_visc_iter",&(E->monitor.max_sdep_visc_iter),"50",m); /* max number of powerlaw iterations */
input_boolean("BDEPV",&(E->viscosity.BDEPV),"off",m);
+
input_int("pdepv_for_flavor",&(E->viscosity.pdepv_for_flavor),"0",m); /* byerlee only for certain flavor */
+ input_int("pdepv_for_zero_comp",&(E->viscosity.pdepv_for_zero_comp),"0",m); /* byerlee only for zero composition */
+
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
input_boolean("cdepv_absolute",&(E->viscosity.cdepv_absolute),"off",m); /* make comp dep viscosity absolute */
@@ -357,7 +361,6 @@
if(E->viscosity.BDEPV)
visc_from_B(E, visc, evisc, propogate);
-
if(E->viscosity.SMOOTH)
apply_viscosity_smoother(E, visc, evisc);
@@ -1258,6 +1261,7 @@
eedot = (float *) safe_malloc((2+nel)*sizeof(float));
+
#ifdef DEBUG
out=fopen("tmp.visc","w");
@@ -1285,7 +1289,10 @@
E->viscosity.plasticity_trans,
E->viscosity.plasticity_viscosity_offset);
fprintf(stderr,"\tpsrw: %i\n",E->viscosity.psrw);
+ if(E->viscosity.pdepv_for_zero_comp && E->viscosity.pdepv_for_flavor)
+ myerror("error, can't have both flavor and composition dependent plasticity",E);
+
if(E->viscosity.pdepv_for_flavor){
if((E->tracers_add_flavors < 1) || (E->control.composition == 0))
myerror("Byerlee is to apply only to certain flavor, but no flavor is set",E);
@@ -1319,20 +1326,24 @@
zz[kk] *= ndz_to_m; /* scale to meters */
if(E->viscosity.pdepv_for_flavor)
tfn[kk]= E->CF[0][node];
+ else if(E->viscosity.pdepv_for_zero_comp)
+ tfn[kk]= E->C[node];
}
for(jj=1;jj <= vpts;jj++) {
zzz = wmax = 0.0;
for(kk=1;kk <= ends;kk++) {
weight = E->N.vpt[GNVINDEX(kk,jj)];
zzz += zz[kk] * weight;
- if(E->viscosity.pdepv_for_flavor){
+ if(E->viscosity.pdepv_for_flavor || E->viscosity.pdepv_for_zero_comp){
if(weight > wmax){
wmax = weight;
point_flavor = tfn[kk]; /* use the most important nodal flavor */
}
}
}
- if((E->viscosity.pdepv_for_flavor==0)||(point_flavor == E->viscosity.pdepv_for_flavor)){
+ if(((E->viscosity.pdepv_for_flavor==0) && (E->viscosity.pdepv_for_zero_comp))||
+ (E->viscosity.pdepv_for_flavor && (point_flavor == E->viscosity.pdepv_for_flavor))||
+ (E->viscosity.pdepv_for_zero_comp && (point_flavor < 0.5))){
npc++;
/* plasticity applies */
if(E->viscosity.plasticity_dimensional){
@@ -1356,7 +1367,9 @@
}
}
}else{
-
+#ifdef DEBUG
+ if(E->parallel.me==0)fprintf(stderr,"calling regular BDEPV %i %i %i %i\n",nel,ends,vpts,visited);
+#endif
/* regular plasticity */
npc=0;
@@ -1376,6 +1389,8 @@
zz[kk] *= ndz_to_m; /* scale to meters */
if(E->viscosity.pdepv_for_flavor) /* nodal flavors */
tfn[kk]= E->CF[0][node];
+ else if(E->viscosity.pdepv_for_zero_comp)
+ tfn[kk]= E->C[node];
}
for(jj=1;jj <= vpts;jj++) {
/* loop over nodes in element */
@@ -1386,14 +1401,30 @@
*/
weight = E->N.vpt[GNVINDEX(kk,jj)];
zzz += zz[kk] * weight;
- if(E->viscosity.pdepv_for_flavor){
+ if(E->viscosity.pdepv_for_flavor || E->viscosity.pdepv_for_zero_comp){
if(weight > wmax){
wmax = weight;
point_flavor = tfn[kk]; /* use the most important nodal flavor */
}
}
}
- if((E->viscosity.pdepv_for_flavor==0)||(point_flavor == E->viscosity.pdepv_for_flavor)){
+ /*
+ three different ways to apply plasticity
+
+ */
+ if(((E->viscosity.pdepv_for_flavor == 0) && (E->viscosity.pdepv_for_zero_comp == 0))|| /* regular operation */
+ (E->viscosity.pdepv_for_flavor && (point_flavor == E->viscosity.pdepv_for_flavor))|| /* flavor
+ dependent
+ plasticity */
+ (E->viscosity.pdepv_for_zero_comp && (point_flavor < 0.5)) /* dependent
+ on
+ regular
+ composition */
+
+ ){
+
+
+ /* */
npc++;
/*
apply plasticity for this integration point if
@@ -1451,7 +1482,7 @@
*/
if(visited)
- fprintf(out,"%10.2f %17.4e %17.4e %17.4e %17.4e %17.4e\n",
+ fprintf(stderr,"%10.2f %17.4e %17.4e %17.4e %17.4e %17.4e\n",
zzz/1e3,tau*E->monitor.tau_scale/1e6,
eedot[i]/E->monitor.time_scale,
ettby*E->data.ref_viscosity,
Modified: mc/3D/CitcomCU/trunk/src/global_defs.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/global_defs.h 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/global_defs.h 2012-04-10 23:55:52 UTC (rev 19936)
@@ -578,6 +578,7 @@
int periodic_x;
int periodic_y;
+ int periodic_pin_or_filter; /* */
double volume;
float layer[4]; /* dimensionless dimensions */
float lidz;
@@ -781,6 +782,7 @@
int adi_heating, visc_heating;
int composition;
+ int composition_neutralize_buoyancy;
float z_comp;
int dfact;
Modified: mc/3D/CitcomCU/trunk/src/prototypes.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/prototypes.h 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/prototypes.h 2012-04-10 23:55:52 UTC (rev 19936)
@@ -397,3 +397,4 @@
int layers(struct All_variables *, float);
int weak_zones(struct All_variables *, int, float);
float boundary_thickness(struct All_variables *, float *);
+void remove_net_drift(struct All_variables *);
Modified: mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2012-04-10 23:51:38 UTC (rev 19935)
+++ mc/3D/CitcomCU/trunk/src/viscosity_descriptions.h 2012-04-10 23:55:52 UTC (rev 19936)
@@ -146,7 +146,7 @@
float pre_comp[CITCOM_CU_VISC_MAXLAYER*2]; /* prefactors */
int layer_pre_comp;
- int pdepv_for_flavor;
+ int pdepv_for_flavor,pdepv_for_zero_comp;
More information about the CIG-COMMITS
mailing list