[cig-commits] r7874 - in mc/3D/CitcomS/trunk: CitcomS/Components/Advection_diffusion lib module

tan2 at geodynamics.org tan2 at geodynamics.org
Thu Aug 23 13:16:10 PDT 2007


Author: tan2
Date: 2007-08-23 13:16:10 -0700 (Thu, 23 Aug 2007)
New Revision: 7874

Modified:
   mc/3D/CitcomS/trunk/CitcomS/Components/Advection_diffusion/Advection_diffusion.py
   mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
   mc/3D/CitcomS/trunk/lib/advection.h
   mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Added new input parameter solver.tsolver.monitor_max_T

When it is on (default), if the max temperature changes too much between timestep,
the temperature field is restored at the tsolver is called using half of the
timestep size.


Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Advection_diffusion/Advection_diffusion.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Advection_diffusion/Advection_diffusion.py	2007-08-23 20:15:50 UTC (rev 7873)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Advection_diffusion/Advection_diffusion.py	2007-08-23 20:16:10 UTC (rev 7874)
@@ -87,6 +87,7 @@
 
         ADV = prop.bool("ADV", default=True)
         filter_temp = prop.bool("filter_temp", default=True)
+        monitor_max_T = prop.bool("monitor_max_T", default=True)
 
         fixed_timestep = prop.float("fixed_timestep", default=0.0)
         finetunedt = prop.float("finetunedt", default=0.9)

Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c	2007-08-23 20:15:50 UTC (rev 7873)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c	2007-08-23 20:16:10 UTC (rev 7874)
@@ -98,25 +98,29 @@
 
   E->advection.timesteps++;
 
-  for(m=1;m<=E->sphere.caps_per_proc;m++)  {
+  for(m=1;m<=E->sphere.caps_per_proc;m++)
     DTdot[m]= (double *)malloc((E->lmesh.nno+1)*sizeof(double));
-    T1[m]= (double *)malloc((E->lmesh.nno+1)*sizeof(double));
-    Tdot1[m]= (double *)malloc((E->lmesh.nno+1)*sizeof(double));
-  }
 
-  for(m=1;m<=E->sphere.caps_per_proc;m++)
-    for (i=1;i<=E->lmesh.nno;i++)   {
-      T1[m][i] = E->T[m][i];
-      Tdot1[m][i] = E->Tdot[m][i];
-    }
 
-  /* get the max temperature for old T */
-  T_interior1 = Tmaxd(E,E->T);
+  if(E->advection.monitor_max_T) {
+     for(m=1;m<=E->sphere.caps_per_proc;m++)  {
+         T1[m]= (double *)malloc((E->lmesh.nno+1)*sizeof(double));
+         Tdot1[m]= (double *)malloc((E->lmesh.nno+1)*sizeof(double));
+     }
 
+     for(m=1;m<=E->sphere.caps_per_proc;m++)
+         for (i=1;i<=E->lmesh.nno;i++)   {
+             T1[m][i] = E->T[m][i];
+             Tdot1[m][i] = E->Tdot[m][i];
+         }
+
+     /* get the max temperature for old T */
+     T_interior1 = Tmaxd(E,E->T);
+  }
+
   E->advection.dt_reduced = 1.0;
   E->advection.last_sub_iterations = 1;
 
-/*   time1= CPU_time0(); */
 
   do {
     E->advection.timestep *= E->advection.dt_reduced;
@@ -134,28 +138,32 @@
       }
     }
 
-    /* get the max temperature for new T */
-    E->monitor.T_interior = Tmaxd(E,E->T);
+    if(E->advection.monitor_max_T) {
+        /* get the max temperature for new T */
+        E->monitor.T_interior = Tmaxd(E,E->T);
 
-    if (E->monitor.T_interior/T_interior1 > E->monitor.T_maxvaried) {
-      for(m=1;m<=E->sphere.caps_per_proc;m++)
-	for (i=1;i<=E->lmesh.nno;i++)   {
-	  E->T[m][i] = T1[m][i];
-	  E->Tdot[m][i] = Tdot1[m][i];
-	}
-      iredo = 1;
-      E->advection.dt_reduced *= 0.5;
-      E->advection.last_sub_iterations ++;
+        /* if the max temperature changes too much, restore the old
+         * temperature field, calling the temperature solver using
+         * half of the timestep size */
+        if (E->monitor.T_interior/T_interior1 > E->monitor.T_maxvaried) {
+            for(m=1;m<=E->sphere.caps_per_proc;m++)
+                for (i=1;i<=E->lmesh.nno;i++)   {
+                    E->T[m][i] = T1[m][i];
+                    E->Tdot[m][i] = Tdot1[m][i];
+                }
+            iredo = 1;
+            E->advection.dt_reduced *= 0.5;
+            E->advection.last_sub_iterations ++;
+        }
     }
 
   }  while ( iredo==1 && E->advection.last_sub_iterations <= 5);
 
+
+  /* filter temperature to remove over-/under-shoot */
   if(E->control.filter_temperature)
     filter(E);
 
-  /*   time0= CPU_time0()-time1; */
-  /*     if(E->control.verbose) */
-  /*       fprintf(E->fp_out,"time=%f\n",time0); */
 
   E->advection.total_timesteps++;
   E->monitor.elapsed_time += E->advection.timestep;
@@ -165,10 +173,15 @@
 
   for(m=1;m<=E->sphere.caps_per_proc;m++) {
     free((void *) DTdot[m] );
-    free((void *) T1[m] );
-    free((void *) Tdot1[m] );
   }
 
+  if(E->advection.monitor_max_T) {
+      for(m=1;m<=E->sphere.caps_per_proc;m++) {
+          free((void *) T1[m] );
+          free((void *) Tdot1[m] );
+      }
+  }
+
   if(E->control.lith_age) {
       if(E->parallel.me==0) fprintf(stderr,"PG_timestep_solve\n");
       lith_age_conform_tbc(E);
@@ -190,6 +203,7 @@
     int m=E->parallel.me;
 
     input_boolean("ADV",&(E->advection.ADVECTION),"on",m);
+    input_boolean("monitor_max_T",&(E->advection.monitor_max_T),"on",m);
 
     input_int("minstep",&(E->advection.min_timesteps),"1",m);
     input_int("maxstep",&(E->advection.max_timesteps),"1000",m);

Modified: mc/3D/CitcomS/trunk/lib/advection.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/advection.h	2007-08-23 20:15:50 UTC (rev 7873)
+++ mc/3D/CitcomS/trunk/lib/advection.h	2007-08-23 20:16:10 UTC (rev 7874)
@@ -44,8 +44,8 @@
   int temp_iterations;
   int sub_iterations;
   int last_sub_iterations;
+  int monitor_max_T;
 
-
 } advection;
 
 

Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c	2007-08-23 20:15:50 UTC (rev 7873)
+++ mc/3D/CitcomS/trunk/module/setProperties.c	2007-08-23 20:16:10 UTC (rev 7874)
@@ -99,6 +99,7 @@
 
     getIntProperty(properties, "ADV", E->advection.ADVECTION, fp);
     getIntProperty(properties, "filter_temp", E->control.filter_temperature, fp);
+    getIntProperty(properties, "monitor_max_T", E->advection.monitor_max_T, fp);
 
     getFloatProperty(properties, "finetunedt", E->advection.fine_tune_dt, fp);
     getFloatProperty(properties, "fixed_timestep", E->advection.fixed_timestep, fp);



More information about the cig-commits mailing list