[cig-commits] r15653 - mc/3D/CitcomCU/trunk/src
becker at geodynamics.org
becker at geodynamics.org
Fri Sep 4 16:32:56 PDT 2009
Author: becker
Date: 2009-09-04 16:32:56 -0700 (Fri, 04 Sep 2009)
New Revision: 15653
Modified:
mc/3D/CitcomCU/trunk/src/Convection.c
mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
mc/3D/CitcomCU/trunk/src/Output_gzdir.c
Log:
Implemented rudimentary restart facility for gzdir output, several
local mods still not in CIG version of CitcomCU.
Modified: mc/3D/CitcomCU/trunk/src/Convection.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Convection.c 2009-09-04 22:50:12 UTC (rev 15652)
+++ mc/3D/CitcomCU/trunk/src/Convection.c 2009-09-04 23:32:56 UTC (rev 15653)
@@ -208,8 +208,11 @@
/* ===============================
Initialization of fields .....
+ NOTE:
+
there's a different routine for GGRD handling in Ggrd_handling
- which has differnet logic for marker init etc.
+ which has differnet logic for marker init etc. this gets called
+ when USE_GGRD is compiled in
=============================== */
@@ -246,7 +249,8 @@
x1 = E->X[1][node];
z1 = E->X[3][node];
y1 = E->X[2][node];
- E->T[node] = 1 - z1 + E->convection.perturb_mag[p] * sin(M_PI * (1.0 - z1)) * cos(E->convection.perturb_k[p] * M_PI * x1) * ((E->mesh.nsd != 3) ? 1.0 : cos(E->convection.perturb_k[p] * M_PI * y1));
+ E->T[node] = 1 - z1 +
+ E->convection.perturb_mag[p] * sin(M_PI * (1.0 - z1)) * cos(E->convection.perturb_k[p] * M_PI * x1) * ((E->mesh.nsd != 3) ? 1.0 : cos(E->convection.perturb_k[p] * M_PI * y1));
// E->T[node] = 1-z1 + E->convection.perturb_mag[p] *
// sin(M_PI*(1.0-z1))*
@@ -289,7 +293,7 @@
0.01*modified_plgndr_a(ll,mm,x1);
*/
- E->T[node] = beta * (1.0 - 1.0 / z1) + E->convection.perturb_mag[p] * modified_plgndr_a(ll, mm, x1) * cos(mm * y1) * sin(M_PI * (z1 - E->sphere.ri) / (E->sphere.ro - E->sphere.ri));
+ //E->T[node] = beta * (1.0 - 1.0 / z1) + E->convection.perturb_mag[p] * modified_plgndr_a(ll, mm, x1) * cos(mm * y1) * sin(M_PI * (z1 - E->sphere.ri) / (E->sphere.ro - E->sphere.ri));
E->T[node] = beta * (1.0 - 1.0 / z1) + E->convection.perturb_mag[p] * sin(M_PI * (E->sphere.ro - z1) / (E->sphere.ro - E->sphere.ri)) * cos(E->convection.perturb_k[p] * M_PI * (x1 - E->XG1[1]) / (E->XG2[1] - E->XG1[1])) * ((E->mesh.nsd != 3) ? 1.0 : cos(E->convection.perturb_k[p] * M_PI * (y1 - E->XG1[2]) / (E->XG2[2] - E->XG1[2])));
@@ -309,9 +313,13 @@
else if(E->control.restart)
{
+#ifdef USE_GZDIR
+ if(E->control.gzdir)
+ process_restart_tc_gzdir(E);
+ else
+#endif
+ process_restart_tc(E);
- process_restart_tc(E);
-
}
temperatures_conform_bcs(E);
Modified: mc/3D/CitcomCU/trunk/src/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2009-09-04 22:50:12 UTC (rev 15652)
+++ mc/3D/CitcomCU/trunk/src/Ggrd_handling.c 2009-09-04 23:32:56 UTC (rev 15653)
@@ -297,8 +297,12 @@
z1 = E->X[3][node];
y1 = E->X[2][node];
/* linear gradient */
- tz = z1* (E->control.TBCtopval-E->control.TBCbotval) +
- E->control.TBCbotval;
+ if(E->mesh.bottbc){
+ tz = z1* (E->control.TBCtopval - E->control.TBCbotval) +
+ E->control.TBCbotval;
+ }else{
+ tz = z1* (E->control.TBCtopval - 1) + 1;
+ }
E->T[node] += tz;
if(E->convection.random_t_init){
@@ -306,7 +310,7 @@
E->T[node] += (-0.5 + drand48()) * E->convection.perturb_mag[0];
}else{
- for(p=0;p<E->convection.number_of_perturbations;p++){ /* perturbations */
+ for(p=0;p < E->convection.number_of_perturbations;p++){ /* perturbations */
E->T[node] += E->convection.perturb_mag[p] *
sin(M_PI * (1.0 - z1)) * cos(E->convection.perturb_k[p] * M_PI * x1) *
@@ -466,8 +470,13 @@
else if(E->control.restart)
{
+#ifdef USE_GZDIR
/* restart part */
- process_restart_tc(E);
+ if(E->control.gzdir)
+ process_restart_tc_gzdir(E);
+ else
+#endif
+ process_restart_tc(E);
}
@@ -482,3 +491,4 @@
return;
}
+
Modified: mc/3D/CitcomCU/trunk/src/Output_gzdir.c
===================================================================
--- mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2009-09-04 22:50:12 UTC (rev 15652)
+++ mc/3D/CitcomCU/trunk/src/Output_gzdir.c 2009-09-04 23:32:56 UTC (rev 15653)
@@ -555,3 +555,81 @@
}
return ((gzFile *)tmp);
}
+
+/* gzdir version, will not limit temperatures to [0;1] */
+void process_restart_tc_gzdir(struct All_variables *E)
+{
+
+ int node, i, j, k, p;
+ FILE *fp;
+ float temp1, temp2, *temp;
+ char input_s[200], input_file[255];
+ gzFile gzin;
+
+ temp = (float *)malloc((E->mesh.noz + 1) * sizeof(float));
+
+ if(E->control.restart == 1 || E->control.restart == 3)
+ {
+
+ sprintf(input_file,"%s/%d/t.%d.%d.gz",
+ E->convection.old_T_file,
+ E->monitor.solution_cycles,E->parallel.me, E->monitor.solution_cycles);
+ gzin = safe_gzopen(input_file, "r");
+ if(gzgets (gzin,input_s, 200) == Z_NULL){
+ fprintf(stderr,"read error\n");
+ parallel_process_termination();
+ }
+ sscanf(input_s, "%i %i %g", &i, &j, &E->monitor.elapsed_time);
+
+ for(node = 1; node <= E->lmesh.nno; node++)
+ {
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%g", &E->T[node]);
+ E->C[node] = 0;
+ E->node[node] = E->node[node] | (INTX | INTZ | INTY);
+ }
+ gzclose(gzin);
+ if(E->parallel.me == 0)
+ fprintf(stderr,"restarting from %s timestep %i time %g\n",
+ E->convection.old_T_file,E->monitor.solution_cycles,
+ E->monitor.elapsed_time);
+
+ }
+ else
+ {
+ fprintf(stderr,"error: restart mode -1 or 2 not implemented for gzdir/vtkout\n");
+ parallel_process_termination();
+ }
+
+ // for composition
+
+ if(E->control.composition && (E->control.restart == 1 || E->control.restart == 2))
+ {
+ convection_initial_markers(E,0);
+ }
+
+ else if(E->control.composition && (E->control.restart == 3))
+ {
+ sprintf(input_file,"%s/%d/c.%d.%d.gz",
+ E->convection.old_T_file, E->monitor.solution_cycles,
+ E->parallel.me, E->monitor.solution_cycles);
+ gzin = safe_gzopen(input_file, "r");
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%d %d ", &i, &j);
+
+ for(node = 1; node <= E->lmesh.nno; node++)
+ {
+ gzgets (gzin,input_s, 200);
+ sscanf(input_s, "%lf", &E->C[node]);
+ }
+ gzclose(gzin);
+
+ convection_initial_markers1(E);
+ }
+
+ E->advection.timesteps = E->monitor.solution_cycles;
+
+ return;
+}
+
+
More information about the CIG-COMMITS
mailing list