[cig-commits] commit: Update the corrections by multiplying or dividing by 1+theta_correction

[Mercurial]

changeset:   59:9d29fb7ccbb6
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Feb 22 05:23:50 2012 -0800
files:       compute_corrections.cxx
description:
Update the corrections by multiplying or dividing by 1+theta_correction


diff -r 32ba145fcf28 -r 9d29fb7ccbb6 compute_corrections.cxx
--- a/compute_corrections.cxx	Wed Feb 22 05:22:31 2012 -0800
+++ b/compute_corrections.cxx	Wed Feb 22 05:23:50 2012 -0800
@@ -1,5 +1,6 @@
 #include "constants.hxx"
 #include "Model.hxx"
+#include <iostream>
 
 extern void compute_v_on_interface(const double zx[Nx+1][Ny],
                                    const double zy[Nx][Ny+1],
@@ -10,6 +11,37 @@ extern void compute_v_on_interface(const
                                    double &vx, double &vy,
                                    double &dvx_y, double &dvy_y, 
                                    double &dvx_yy, double &dvy_yy);
+
+
+void update_correction(const double &current, const double &theta,
+                       double &old)
+{
+  if(old==0)
+    old=current;
+  else
+    {
+      if(current/old>1+theta)
+        old*=1+theta;
+      else if(current/old<1-theta)
+        {
+          if(current*old>1)
+            old*=1-theta;
+          else
+            {
+              double temp(old-(old-current)*theta);
+              if(temp*old<0)
+                {
+                  std::cout << "switch\n";
+                old=0;
+                }
+              else
+                old=temp;
+            }
+        }
+      else
+        old=current;
+    }
+}
 
 void compute_corrections(const Model &model, const double zx[Nx+1][Ny],
                          const double zy[Nx][Ny+1],
@@ -62,7 +94,9 @@ void compute_corrections(const Model &mo
                     Cx_new+=(p_jump + dx*dp_x_jump)/h;
                     Cx_new-=eta_jump*(dvy_y - dx*dvx_yy)/h;
                   }
-                Cx[i][j]+=theta_correction*(Cx_new-Cx[i][j]);
+
+                update_correction(Cx_new,theta_correction,Cx[i][j]);
+                // Cx[i][j]+=theta_correction*(Cx_new-Cx[i][j]);
               }
 
             /* Cy */
@@ -97,7 +131,8 @@ void compute_corrections(const Model &mo
                       eta_jump*(dvx_y - dx*dvy_yy)/h;
                     Cy_new-=dzx_yx_correction;
                   }
-                Cy[i][j]+=theta_correction*(Cy_new-Cy[i][j]);
+                update_correction(Cy_new,theta_correction,Cy[i][j]);
+                // Cy[i][j]+=theta_correction*(Cy_new-Cy[i][j]);
               }
 
             /* Cp */
@@ -118,7 +153,8 @@ void compute_corrections(const Model &mo
                 double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
                 Cp_new=-dzx_x_correction;
 
-                Cp[i][j]+=theta_correction*(Cp_new-Cp[i][j]);
+                update_correction(Cp_new,theta_correction,Cp[i][j]);
+                // Cp[i][j]+=theta_correction*(Cp_new-Cp[i][j]);
               }
           }
     }