[cig-commits] commit: Fix a factor of 2 bug in the second derivative when adding an element to the interface interpolation matrix

[Mercurial]

changeset:   47:06750548dbd3
user:        Walter Landry <wlandry at caltech.edu>
date:        Tue Feb 14 16:23:18 2012 -0800
files:       compute_v_on_interface.cxx
description:
Fix a factor of 2 bug in the second derivative when adding an element to the interface interpolation matrix


diff -r fdae064362c2 -r 06750548dbd3 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx	Mon Feb 13 20:06:14 2012 -0800
+++ b/compute_v_on_interface.cxx	Tue Feb 14 16:23:18 2012 -0800
@@ -1,6 +1,7 @@
 #include <algorithm>
 #include <cmath>
 #include <iostream>
+#include <iomanip>
 #include "constants.hxx"
 #include <map>
 #include <gsl/gsl_linalg.h>
@@ -24,8 +25,18 @@ namespace {
     m[5+n*n_terms]+=sign*dxm*dxm/2;
     m[6+n*n_terms]+=sign*dxp;
     m[7+n*n_terms]+=sign*dxp*dy;
-    m[8+n*n_terms]+=sign*dxp*dxp;
+    m[8+n*n_terms]+=sign*dxp*dxp/2;
   }
+}
+
+double analytic(const double &x, const double &y)
+{
+  double result=1.1 + 1.2*y + 1.3*y*y/2 + 1.4*y*y*y/6;
+  if(x>0)
+    return result
+      + 1.5*x + 1.6*x*x/2 + 1.7*x*x*x/6 + 1.8*x*y + 1.9*x*x*y/2 + 2.0*x*y*y/2;
+  return result
+    + 2.5*x + 2.6*x*x/2 + 2.7*x*x*x/6 + 2.8*x*y + 2.9*x*x*y/2 + 3.0*x*y*y/2;
 }
 
 template<int Ny>
@@ -72,22 +83,45 @@ void interpolate_to_surface(double z[][N
          i_center+1,j+1 - i_center-1,j-1
       */
 
-      add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,0,-di*dx,0,m,rhs);
-      add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,0,-di*dx,h,m,rhs);
-      add_to_matrix(n_terms,2,z,log_eta,i_center,j-1,1,0,-di*dx,-h,m,rhs);
-      add_to_matrix(n_terms,3,z,log_eta,i_center-di,j,1,0,-di*(dx+h),0,m,rhs);
-      add_to_matrix(n_terms,4,z,log_eta,i_center-2*di,j,1,0,-di*(dx+2*h),
+      add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,0,di*dx,0,m,rhs);
+      add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,0,di*dx,h,m,rhs);
+      add_to_matrix(n_terms,2,z,log_eta,i_center,j-1,1,0,di*dx,-h,m,rhs);
+      add_to_matrix(n_terms,3,z,log_eta,i_center-di,j,1,0,di*(dx+h),0,m,rhs);
+      add_to_matrix(n_terms,4,z,log_eta,i_center-2*di,j,1,0,di*(dx+2*h),
                     0,m,rhs);
-      add_to_matrix(n_terms,5,z,log_eta,i_center-di,j+1,1,0,-di*(dx+h),h,m,rhs);
-      add_to_matrix(n_terms,5,z,log_eta,i_center-di,j-1,1,0,-di*(dx+h),
+      add_to_matrix(n_terms,5,z,log_eta,i_center-di,j+1,1,0,di*(dx+h),h,m,rhs);
+      add_to_matrix(n_terms,5,z,log_eta,i_center-di,j-1,1,0,di*(dx+h),
                     -h,m,rhs,-1);
 
-      add_to_matrix(n_terms,6,z,log_eta,i_center+di,j,1,-di*(dx-h),0,0,m,rhs);
-      add_to_matrix(n_terms,7,z,log_eta,i_center+2*di,j,1,-di*(dx-2*h),0,
+      add_to_matrix(n_terms,6,z,log_eta,i_center+di,j,1,di*(dx-h),0,0,m,rhs);
+      add_to_matrix(n_terms,7,z,log_eta,i_center+2*di,j,1,di*(dx-2*h),0,
                     0,m,rhs);
-      add_to_matrix(n_terms,8,z,log_eta,i_center+di,j+1,1,-di*(dx-h),0,h,m,rhs);
-      add_to_matrix(n_terms,8,z,log_eta,i_center+di,j-1,1,-di*(dx-h),0,
+      add_to_matrix(n_terms,8,z,log_eta,i_center+di,j+1,1,di*(dx-h),0,h,m,rhs);
+      add_to_matrix(n_terms,8,z,log_eta,i_center+di,j-1,1,di*(dx-h),0,
                     -h,m,rhs,-1);
+
+
+      // if(j>0)
+      //   {
+      // std::cout << "Set rhs "
+      //           << di << " "
+      //           << dx << " "
+      //           << h << " "
+      //           << dx+h << " "
+      //           << dx-h << " "
+      //           << dx-2*h << " "
+      //           << "\n";
+      // rhs[0]=analytic(di*dx,0);
+      // rhs[1]=analytic(di*dx,h);
+      // rhs[2]=analytic(di*dx,-h);
+      // rhs[3]=analytic(di*(dx+h),0);
+      // rhs[4]=analytic(di*(dx+2*h),0);
+      // rhs[5]=analytic(di*(dx+h),h)-analytic(di*(dx+h),-h);
+
+      // rhs[6]=analytic(di*(dx-h),0);
+      // rhs[7]=analytic(di*(dx-2*h),0);
+      // rhs[8]=analytic(di*(dx-h),h)-analytic(di*(dx-h),-h);
+      //   }
 
       /* At the j==0 boundary, all of the j-1 terms are not applied.
          Instead, for n=2, apply dvx/dy=0. */
@@ -139,13 +173,17 @@ void interpolate_to_surface(double z[][N
                     h/2,m,rhs);
     }
 
-  // if(j>510)
+  // if(j>0)
   //   for(int a=0;a<9;++a)
   //     {
   //       std::cout << a << " "
+  //                 << std::setw(15)
+  //                 << std::setprecision(10)
   //                 << rhs[a] << " ";
   //       for(int b=0;b<9;++b)
-  //         std::cout << m[b+9*a] << " ";
+  //         std::cout << std::setw(15)
+  //                   << std::setprecision(10)
+  //                   << m[b+9*a] << " ";
   //       std::cout << "\n";
   //     }
 
@@ -156,6 +194,15 @@ void interpolate_to_surface(double z[][N
   gsl_permutation *perm=gsl_permutation_alloc(n_terms);
   gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
   gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,result);
+
+  // if(j>0)
+  //   {
+  // std::cout << "vector ";
+  // for(int ii=0;ii<9;++ii)
+  //   std::cout << gsl_vector_get(result,ii) << " ";
+  // std::cout << "\n";
+  // exit(0);
+  //   }
 
   v0=gsl_vector_get(result,0);
   dv_y=gsl_vector_get(result,1);