[cig-commits] commit: Split out interface initialization into compute_v_on_interface
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:26 PST 2012
changeset: 27:bdb505cb7635
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 07 03:28:44 2012 -0800
files: compute_v_on_interface.cxx main.cxx wscript
description:
Split out interface initialization into compute_v_on_interface
diff -r 99f0c4630c8c -r bdb505cb7635 compute_v_on_interface.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface.cxx Tue Feb 07 03:28:44 2012 -0800
@@ -0,0 +1,96 @@
+#include <gsl/gsl_spline.h>
+#include "constants.hxx"
+
+template<class T>
+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]);
+ double um=z[i-1][j]*exp(-log_eta[i-1][j]);
+
+ double du=(up-um)/(2*h);
+ double ddu=(up - 2*u0 + um)/(h*h);
+
+ return u0 + x*du + x*x*ddu/2;
+}
+
+void compute_v_on_interface(double zx[N+1][N], double zy[N][N+1],
+ double log_etax[N+1][N], double log_etay[N][N+1],
+ double vx[N], double dvx_y[N+1], double dvx_yy[N],
+ double vy[N+1], double dvy_y[N], double dvy_yy[N+1])
+{
+ 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;
+
+ /* vx, dvx_y, dvx_yy */
+ for(int j=0;j<N;++j)
+ {
+ double low1=(zx[i][j]*exp(-log_etax[i][j])*(dx+1)
+ - zx[i-1][j]*exp(-log_etax[i-1][j])*dx);
+
+ 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;
+ }
+
+ 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);
+
+ 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);
+ double dx=middle/h-0.5-i;
+
+ /* vy, dvy_y, dvy_yy */
+ for(int j=0;j<N+1;++j)
+ {
+ double low1=(zy[i][j]*exp(-log_etay[i][j])*(dx+1)
+ - zy[i-1][j]*exp(-log_etay[i-1][j])*dx);
+
+ double high1=(zy[i+1][j]*exp(-log_etay[i+1][j])*(2-dx)
+ - zy[i+2][j]*exp(-log_etay[i+2][j])*(1-dx));
+
+ 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;
+ }
+ 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);
+
+ 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);
+
+ for(int i=0;i<N;++i)
+ dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
+
+ gsl_spline_free(vy_spline);
+ gsl_interp_accel_free(vy_accel);
+ }
+}
diff -r 99f0c4630c8c -r bdb505cb7635 main.cxx
--- a/main.cxx Tue Feb 07 03:10:41 2012 -0800
+++ b/main.cxx Tue Feb 07 03:28:44 2012 -0800
@@ -1,7 +1,6 @@
#include <iostream>
#include <fstream>
#include <sstream>
-#include <gsl/gsl_spline.h>
#include "Model.hxx"
#include "constants.hxx"
@@ -11,42 +10,13 @@ extern void initial(const Model &model,
double log_etax[N+1][N], double log_etay[N][N+1],
double p[N][N], double fx[N+1][N], double fy[N][N+1]);
-
-template<class T>
-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]);
- double um=z[i-1][j]*exp(-log_eta[i-1][j]);
-
- double du=(up-um)/(2*h);
- double ddu=(up - 2*u0 + um)/(h*h);
-
- 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);
-// }
+extern void compute_v_on_interface(double zx[N+1][N], double zy[N][N+1],
+ double log_etax[N+1][N],
+ double log_etay[N][N+1],
+ double vx[N], double dvx_y[N+1],
+ double dvx_yy[N],
+ double vy[N+1], double dvy_y[N],
+ double dvy_yy[N+1]);
int main()
{
@@ -85,79 +55,8 @@ 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;
-
- /* vx, dvx_y, dvx_yy */
- for(int j=0;j<N;++j)
- {
- double low1=(zx[i][j]*exp(-log_etax[i][j])*(dx+1)
- - zx[i-1][j]*exp(-log_etax[i-1][j])*dx);
-
- 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;
- }
-
- 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);
-
- 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);
- double dx=middle/h-0.5-i;
-
- /* vy, dvy_y, dvy_yy */
- for(int j=0;j<N+1;++j)
- {
- double low1=(zy[i][j]*exp(-log_etay[i][j])*(dx+1)
- - zy[i-1][j]*exp(-log_etay[i-1][j])*dx);
-
- double high1=(zy[i+1][j]*exp(-log_etay[i+1][j])*(2-dx)
- - zy[i+2][j]*exp(-log_etay[i+2][j])*(1-dx));
-
- 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;
- }
- 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);
-
- 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);
-
- for(int i=0;i<N;++i)
- dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
-
- gsl_spline_free(vy_spline);
- gsl_interp_accel_free(vy_accel);
- }
+ compute_v_on_interface(zx,zy,log_etax,log_etay,vx,dvx_y,dvx_yy,
+ vy,dvy_y,dvy_yy);
}
/* Smoothing sweeps */
diff -r 99f0c4630c8c -r bdb505cb7635 wscript
--- a/wscript Tue Feb 07 03:10:41 2012 -0800
+++ b/wscript Tue Feb 07 03:28:44 2012 -0800
@@ -10,7 +10,8 @@ def build(bld):
bld.program(
features = ['cxx','cprogram'],
source = ['main.cxx',
- 'initial.cxx'],
+ 'initial.cxx',
+ 'compute_v_on_interface.cxx'],
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG'],
More information about the CIG-COMMITS
mailing list