[cig-commits] commit: Use splines to compute derivatives

[Mercurial]

changeset:   26:99f0c4630c8c
user:        Walter Landry <wlandry at caltech.edu>
date:        Tue Feb 07 03:10:41 2012 -0800
files:       constants.hxx main.cxx
description:
Use splines to compute derivatives


diff -r b28d5b89ebe7 -r 99f0c4630c8c constants.hxx
--- a/constants.hxx	Mon Feb 06 15:33:20 2012 -0800
+++ b/constants.hxx	Tue Feb 07 03:10:41 2012 -0800
@@ -1,12 +1,12 @@ const int N(1024);
 const int N(1024);
 const double min_eta=1;
-const double max_eta=2;
+const double max_eta=1e4;
 const double eta_jump=max_eta-min_eta;
 #include <cmath>
 const double log_max_eta=std::log(max_eta);
 const double log_min_eta=std::log(min_eta);
-// const double middle=25.0/64 + 0.0000000001;
-const double middle=0.4;
+const double middle=25.00000001/64;
+// const double middle=0.4;
 const double h(1.0/N);
 const double pi(atan(1.0)*4);
 
diff -r b28d5b89ebe7 -r 99f0c4630c8c main.cxx
--- a/main.cxx	Mon Feb 06 15:33:20 2012 -0800
+++ b/main.cxx	Tue Feb 07 03:10:41 2012 -0800
@@ -1,6 +1,7 @@
 #include <iostream>
 #include <fstream>
 #include <sstream>
+#include <gsl/gsl_spline.h>
 
 #include "Model.hxx"
 #include "constants.hxx"
@@ -12,8 +13,8 @@ extern void initial(const Model &model, 
 
 
 template<class T>
-double interpolate_v(T z, T log_eta,
-                     const int &i, const int &j, const double &h, const double &x)
+double interpolate_v(T z, T log_eta, const int &i, const int &j,
+                     const double &h, const double &x)
 {
   double u0=z[i][j]*exp(-log_eta[i][j]);
   double up=z[i+1][j]*exp(-log_eta[i+1][j]);
@@ -24,6 +25,28 @@ double interpolate_v(T z, T log_eta,
 
   return u0 + x*du + x*x*ddu/2;
 }
+
+// template<class T>
+// void deriv(T v, T dv, const int n)
+// {
+//   dv[0]=dv[n]=0;
+//   dv[1]=(v[1]-v[0])/h;
+//   dv[n-1]=(v[n-1]-v[n-2])/h;
+
+//   for(int j=2;j<n-1;++j)
+//     dv[j]=((v[j]-v[j-1])*9/8 - (v[j+1]-v[j-2])/24)/h;
+// }
+
+// template<class T>
+// void deriv_2(T v, T ddv, const int n)
+// {
+//   ddv[0]=ddv[n]=0;
+//   ddv[1]=(v[2]-2*v[1]+v[0])/(h*h);
+//   ddv[n-1]=(v[n]-2*v[n-1]+v[n-2])/(h*h);
+
+//   for(int j=2;j<n-1;++j)
+//     ddv[j]=(-3.5*v[j] + 2*(v[j+1]+v[j-1]) -.25*(v[j+2]+v[j-2]))/(h*h);
+// }
 
 int main()
 {
@@ -62,6 +85,12 @@ int main()
 
       if(model==Model::solCx)
         {
+          double y_x[N], y_y[N+1];
+          for(int i=0;i<N;++i)
+            y_x[i]=h*(i+0.5);
+          for(int i=0;i<N+1;++i)
+            y_y[i]=h*i;
+
           {
             int i(middle/h);
             double dx=middle/h-i;
@@ -75,19 +104,26 @@ int main()
                 double high1=(zx[i+1][j]*exp(-log_etax[i+1][j])*(2-dx)
                              - zx[i+2][j]*exp(-log_etax[i+2][j])*(1-dx));
 
-
                 double low=interpolate_v(zx,log_etax,i-1,j,h,(1+dx)*h);
                 double high=interpolate_v(zx,log_etax,i+2,j,h,(-2+dx)*h);
 
                 vx[j]=low*(1-dx) + high*dx;
               }
 
-            dvx_y[0]=dvx_y[N]=0;
-            for(int j=1;j<N;++j)
-              dvx_y[j]=(vx[j]-vx[j-1])/h;
+            gsl_interp_accel *vx_accel=gsl_interp_accel_alloc ();
+            gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N);
+            gsl_spline_init(vx_spline, y_x, vx, N);
 
-            for(int j=0;j<N;++j)
-              dvx_yy[j]=(dvx_y[j+1]-dvx_y[j])/h;
+            dvx_y[0]=0;
+            dvx_y[N]=0;
+            for(int i=1;i<N;++i)
+              dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
+
+            for(int i=0;i<N;++i)
+              dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i],vx_accel);
+
+            gsl_spline_free(vx_spline);
+            gsl_interp_accel_free(vx_accel);            
           }
           {
             int i(middle/h-0.5);
@@ -105,36 +141,22 @@ int main()
                 double low=interpolate_v(zy,log_etay,i-1,j,h,h*(dx+1));
                 double high=interpolate_v(zy,log_etay,i+2,j,h,h*(-2+dx));
 
-                // vy[j]=low*(1-dx) + high*dx;
+                vy[j]=low*(1-dx) + high*dx;
+              }
+            gsl_interp_accel *vy_accel=gsl_interp_accel_alloc ();
+            gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline, N+1);
+            gsl_spline_init(vy_spline, y_y, vy, N+1);
 
-                vy[j]=high;
+            dvy_y[0]=0;
+            dvy_y[N]=0;
+            for(int i=1;i<N+1;++i)
+              dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i],vy_accel);
 
-              if(j>=N/4-2 && j<=N/4+2)
-                // if(j==N/4)
-                  std::cout << "vy "
-                            << j << " "
-                            << vy[j] << " "
-                            << low << " "
-                            << high << " "
-                            << log_etay[i][j] << " "
-                            << log_etay[i+1][j] << " "
-                            << log_etay[i-1][j] << " "
-                            // << dx << " "
-                            // << zy[i-2][j]*exp(-log_etay[i-2][j]) << " "
-                            // << zy[i-1][j]*exp(-log_etay[i-1][j]) << " "
-                            // << zy[i][j]*exp(-log_etay[i][j]) << " "
-                            // << zy[i+1][j]*exp(-log_etay[i+1][j]) << " "
-                            // << zy[i+2][j]*exp(-log_etay[i+2][j]) << " "
-                            // << zy[i+3][j]*exp(-log_etay[i+3][j]) << " "
-                            << "\n";
-              }
+            for(int i=0;i<N;++i)
+              dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
 
-            for(int j=0;j<N;++j)
-              dvy_y[j]=(vy[j+1]-vy[j])/h;
-
-            dvy_yy[0]=dvy_yy[N]=0;
-            for(int j=1;j<N;++j)
-              dvy_yy[j]=(dvy_y[j+1]-dvy_y[j])/h;
+            gsl_spline_free(vy_spline);
+            gsl_interp_accel_free(vy_accel);            
           }
         }