[cig-commits] commit: Use Nx and Ny instead of just N
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:37 PST 2012
changeset: 36:127bfe48757b
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 08 03:27:21 2012 -0800
files: compute_v_on_interface.cxx constants.hxx initial.cxx main.cxx
description:
Use Nx and Ny instead of just N
diff -r d818fb4febda -r 127bfe48757b compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Wed Feb 08 01:35:41 2012 -0800
+++ b/compute_v_on_interface.cxx Wed Feb 08 03:27:21 2012 -0800
@@ -19,27 +19,29 @@ double interpolate_v(T z, T log_eta, con
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+4], double dvx_y[N+1], double dvx_yy[N],
- double vy[N+1], double dvy_y[N], double dvy_yy[N+1])
+void compute_v_on_interface(double zx[Nx+1][Ny], double zy[Nx][Ny+1],
+ double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
+ double vx[Ny+4], double dvx_y[Ny+1],
+ double dvx_yy[Ny],
+ double vy[Ny+1], double dvy_y[Ny],
+ double dvy_yy[Ny+1])
{
- double y_x[N+4], y_y[N+1];
- for(int i=0;i<N+4;++i)
+ double y_x[Ny+4], y_y[Ny+1];
+ for(int i=0;i<Ny+4;++i)
y_x[i]=h*(i-1.5);
- for(int i=0;i<N+1;++i)
+ for(int i=0;i<Ny+1;++i)
y_y[i]=h*i;
const double a1(h), a2(2*a1), b1(a1), b2(a2);
- double ux1_p[N], ux2_p[N], ux1_m[N], ux2_m[N],
- uy1_p[N+1], uy2_p[N+1], uy1_m[N+1], uy2_m[N+1];
+ double ux1_p[Ny], ux2_p[Ny], ux1_m[Ny], ux2_m[Ny],
+ uy1_p[Ny+1], uy2_p[Ny+1], uy1_m[Ny+1], uy2_m[Ny+1];
{
int i(middle/h);
double dx=middle/h-i;
/* vx, dvx_y, dvx_yy */
- for(int j=0;j<N;++j)
+ for(int j=0;j<Ny;++j)
{
double low1=(zx[i][j]*exp(-log_etax[i][j])*(dx+1)
- zx[i-1][j]*exp(-log_etax[i-1][j])*dx);
@@ -63,21 +65,21 @@ void compute_v_on_interface(double zx[N+
+ zx[i-1][j]*exp(-log_etax[i-1][j])*dx;
}
// vx[0]=vx[3];
- // vx[N+3]=vx[N];
+ // vx[Ny+3]=vx[Ny];
// vx[1]=vx[2];
- // vx[N+2]=vx[N+1];
+ // vx[Ny+2]=vx[Ny+1];
- vx[0]=-vx[N];
- vx[N+3]=-vx[3];
- vx[1]=-vx[N+1];
- vx[N+2]=-vx[2];
+ vx[0]=-vx[Ny];
+ vx[Ny+3]=-vx[3];
+ vx[1]=-vx[Ny+1];
+ vx[Ny+2]=-vx[2];
}
{
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)
+ for(int j=0;j<Ny+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);
@@ -114,15 +116,15 @@ void compute_v_on_interface(double zx[N+
/* Compute the derivatives from the current guess */
/* dvx_dy, dvx_dyy */
gsl_interp_accel *vx_accel=gsl_interp_accel_alloc ();
- gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, N+4);
- gsl_spline_init(vx_spline, y_x, vx, N+4);
+ gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline, Ny+4);
+ gsl_spline_init(vx_spline, y_x, vx, Ny+4);
dvx_y[0]=0;
- dvx_y[N]=0;
- for(int i=1;i<N;++i)
+ dvx_y[Ny]=0;
+ for(int i=1;i<Ny;++i)
dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
- for(int i=0;i<N;++i)
+ for(int i=0;i<Ny;++i)
dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i+2],vx_accel);
gsl_spline_free(vx_spline);
@@ -130,20 +132,24 @@ void compute_v_on_interface(double zx[N+
/* dvy_dy, dvy_dyy */
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);
+ gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline, Ny+1);
+ gsl_spline_init(vy_spline, y_y, vy, Ny+1);
- for(int i=0;i<N;++i)
+ for(int i=0;i<Ny;++i)
dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+2],vy_accel);
- for(int i=0;i<N;++i)
+ for(int i=0;i<Ny;++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);
+ // std::cout << "vx "
+ // << vx[Ny/4+2] << " "
+ // << "\n";
+
/* Update the current guess based on the derivatives */
- for(int i=0;i<N;++i)
+ for(int i=0;i<Ny;++i)
{
double zx_jump=-eta_jump*dvy_y[i];
double temp=
@@ -159,16 +165,16 @@ void compute_v_on_interface(double zx[N+
vx[i+2]=(1-relax)*vx[i+2] + relax*temp;
}
// vx[0]=vx[3];
- // vx[N+3]=vx[N];
+ // vx[Ny+3]=vx[Ny];
// vx[1]=vx[2];
- // vx[N+2]=vx[N+1];
+ // vx[Ny+2]=vx[Ny+1];
- vx[0]=-vx[N];
- vx[N+3]=-vx[3];
- vx[1]=-vx[N+1];
- vx[N+2]=-vx[2];
+ vx[0]=-vx[Ny];
+ vx[Ny+3]=-vx[3];
+ vx[1]=-vx[Ny+1];
+ vx[Ny+2]=-vx[2];
- for(int i=0;i<N+1;++i)
+ for(int i=0;i<Ny+1;++i)
{
double zy_jump=-eta_jump*dvx_y[i];
double temp=
diff -r d818fb4febda -r 127bfe48757b constants.hxx
--- a/constants.hxx Wed Feb 08 01:35:41 2012 -0800
+++ b/constants.hxx Wed Feb 08 03:27:21 2012 -0800
@@ -1,4 +1,6 @@ const int N(128);
-const int N(128);
+const int NN(1024);
+const int Nx(NN);
+const int Ny(NN);
const double min_eta=1;
const double max_eta=1e6;
const double eta_jump=max_eta-min_eta;
@@ -7,7 +9,7 @@ const double log_min_eta=std::log(min_et
const double log_min_eta=std::log(min_eta);
const double middle=25.00000001/64;
// const double middle=0.4;
-const double h(1.0/N);
+const double h(1.0/Nx);
const double pi(atan(1.0)*4);
const double r2_inclusion=0.1;
diff -r d818fb4febda -r 127bfe48757b initial.cxx
--- a/initial.cxx Wed Feb 08 01:35:41 2012 -0800
+++ b/initial.cxx Wed Feb 08 03:27:21 2012 -0800
@@ -6,9 +6,9 @@
#include <cmath>
#include <gsl/gsl_linalg.h>
-void initial(const Model &model, double zx[N+1][N], double zy[N][N+1],
- 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])
+void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
+ double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1])
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
@@ -36,17 +36,6 @@ void initial(const Model &model, double
(kx0*sinh0 + kx1*sinh1)*rho/2,
((sinh0 + kx0*cosh0) + (sinh1 + kx1*cosh1))*rho/2};
- // double m_temp[]={m[0],m[1],m[2],m[3],m[4],m[5],m[6],m[7],m[8],m[9],m[10],m[11],m[12],m[13],m[14],m[15],rhs[3]};
-
- // std::cout << "rhs "
- // << rhs[0] << " "
- // << rhs[1] << " "
- // << rhs[2] << " "
- // << rhs[3] << " "
- // << (sinh0 + kx0*cosh0)*rho/2 << " "
- // << (sinh1 + kx1*cosh1)*rho/2 << " "
- // << "\n";
-
gsl_matrix_view mv=gsl_matrix_view_array(m,4,4);
gsl_vector_view rhsv=gsl_vector_view_array(rhs,4);
gsl_vector *x=gsl_vector_alloc(4);
@@ -60,17 +49,10 @@ void initial(const Model &model, double
gsl_vector_free(x);
gsl_permutation_free(perm);
- // std::cout << "v tau "
- // << v0 << " "
- // << tau_eta_0 << " "
- // << v1 << " "
- // << tau_eta_1 << " "
- // << "\n";
-
- for(int i=0;i<N+1;++i)
- for(int j=0;j<N+1;++j)
+ for(int i=0;i<Nx+1;++i)
+ for(int j=0;j<Ny+1;++j)
{
- if(i<N && j<N)
+ if(i<Nx && j<Ny)
{
double x((i+0.5)*h), y((j+0.5)*h), r2(x*x+y*y);
switch(model)
@@ -118,7 +100,7 @@ void initial(const Model &model, double
break;
}
}
- if(j<N)
+ if(j<Ny)
{
double x(i*h), y((j+0.5)*h), r2(x*x+y*y);
zx[i][j]=0;
@@ -174,25 +156,6 @@ void initial(const Model &model, double
double zy=(v*csh + (v+tau_eta)*kx*snh
- sign*rho*(cs - csh - kx*snh)/(2*eta))
*std::sin(pi*y)*eta;
-
- // if(j==N/4 && i>middle*N-5 && i<middle*N+5)
- // std::cout << "zx "
- // << i << " "
- // << j << " "
- // << x << " "
- // << y << " "
- // // << tau_xx << " "
- // // << dvx_z << " "
- // // << dvz_x << " "
- // // << dvx_z + dvz_x << " "
- // // << tau_xz << " "
- // // << zx[i][j] << " "
- // << zx[i][j]*2*pi*pi*eta_jump/eta << " "
- // // << 2*zy*pi*pi*eta_jump/eta << " "
- // << eta_jump*2*zy*pi/(tan(pi*y)*eta) << " "
- // // << p[i][j] << " "
- // << zy/eta << " "
- // << "\n";
}
break;
case Model::sinker:
@@ -227,7 +190,7 @@ void initial(const Model &model, double
}
}
- if(i<N)
+ if(i<Nx)
{
double x((i+0.5)*h), y(j*h), r2(x*x+y*y);
zy[i][j]=0;
@@ -272,16 +235,6 @@ void initial(const Model &model, double
zy[i][j]=(v*csh + (v+tau_eta)*kx*snh
- sign*rho*(cs - csh - kx*snh)/(2*eta))
*std::sin(pi*y)*eta;
-
- // if(j==N/4 && i>middle*N-2 && i<middle*N+2)
- // std::cout << "zy "
- // << i << " "
- // << j << " "
- // << x << " "
- // << y << " "
- // << zy[i][j]/eta << " "
- // << "\n";
-
}
break;
case Model::sinker:
@@ -318,9 +271,9 @@ void initial(const Model &model, double
}
}
}
- write_vtk("p_initial",N,N,p);
- write_vtk("zx_initial",N+1,N,zx);
- write_vtk("zy_initial",N,N,zy);
- write_vtk("log_etax",N+1,N,log_etax);
- write_vtk("log_etay",N,N,log_etay);
+ write_vtk("p_initial",Nx,Ny,p);
+ write_vtk("zx_initial",Nx+1,Ny,zx);
+ write_vtk("zy_initial",Nx,Ny,zy);
+ write_vtk("log_etax",Nx+1,Ny,log_etax);
+ write_vtk("log_etay",Nx,Ny,log_etay);
}
diff -r d818fb4febda -r 127bfe48757b main.cxx
--- a/main.cxx Wed Feb 08 01:35:41 2012 -0800
+++ b/main.cxx Wed Feb 08 03:27:21 2012 -0800
@@ -6,17 +6,17 @@
#include "constants.hxx"
#include "write_vtk.hxx"
-extern void initial(const Model &model, double zx[N+1][N], double zy[N][N+1],
- 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]);
+extern void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
+ double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1]);
-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+4], double dvx_y[N+1],
- double dvx_yy[N],
- double vy[N+1], double dvy_y[N],
- double dvy_yy[N+1]);
+extern void compute_v_on_interface(double zx[Nx+1][Ny], double zy[Nx][Ny+1],
+ double log_etax[Nx+1][Ny],
+ double log_etay[Nx][Ny+1],
+ double vx[Ny+4], double dvx_y[Ny+1],
+ double dvx_yy[Ny],
+ double vy[Ny+1], double dvy_y[Ny],
+ double dvy_yy[Ny+1]);
int main()
{
@@ -26,8 +26,8 @@ int main()
the continuity residual more messy. Also, with it collocated
with z, it is easy to compute v=z/eta. */
- double zx[N+1][N], zy[N][N+1], log_etax[N+1][N], log_etay[N][N+1],
- p[N][N], dp[N][N], div[N][N], fx[N+1][N], fy[N][N+1];
+ double zx[Nx+1][Ny], zy[Nx][Ny+1], log_etax[Nx+1][Ny], log_etay[Nx][Ny+1],
+ p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1];
/* Initial conditions */
@@ -41,7 +41,7 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
- double Resid_p[N][N], Resid_x[N+1][N], Resid_y[N][N+1];
+ double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
double dzx_x_correction, dzx_xx_correction, dzx_yx_correction,
dzy_xx_correction, dzy_xy_correction;
@@ -50,8 +50,8 @@ int main()
{
/* Compute the boundary jumps */
- double vx[N+4], dvx_y[N+1], dvx_yy[N],
- vy[N+1], dvy_y[N], dvy_yy[N+1];
+ double vx[Ny+4], dvx_y[Ny+1], dvx_yy[Ny],
+ vy[Ny+1], dvy_y[Ny], dvy_yy[Ny+1];
if(model==Model::solCx)
{
@@ -69,19 +69,20 @@ int main()
/* zx */
for(int rb=0;rb<2;++rb)
{
- for(int j=0;j<N;++j)
+ for(int j=0;j<Ny;++j)
{
int i_min=(j+rb)%2;
- for(int i=i_min;i<N+1;i+=2)
+ for(int i=i_min;i<Nx+1;i+=2)
{
- if(i!=0 && i!=N)
+ if(i!=0 && i!=Nx)
{
/* Do the finite difference parts */
const double x(i*h), y((j+0.5)*h);
/* Derivatives of z */
- double dzx_xx=(zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
+ double dzx_xx=
+ (zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
- (zy[i][j] - zy[i-1][j]))/(h*h);
@@ -96,14 +97,15 @@ int main()
dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
}
- else if(j==N-1)
+ else if(j==Ny-1)
{
dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
}
else
{
- dzx_yy=(zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ dzx_yy=
+ (zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
dzx_y=(zx[i][j+1] - zx[i][j-1])/(2*h);
}
@@ -145,148 +147,57 @@ int main()
dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
dlog_etaxy=0;
- // if(j==N/4)
- // std::cout << "dzx_xx "
- // << i << " "
- // << j << " "
- // << x << " "
- // << y << " "
- // // << zx[i][j] << " "
- // << zy[i][j] << " "
- // // << log_etax[i][j] << " "
- // // << log_etay[i][j] << " "
- // << "\n";
-
if(x-h<middle && x+h>middle)
{
double dx=(x>middle) ? (middle-(x-h)) : (middle-(x+h));
double zx_jump=eta_jump*vx[j+2];
double dzx_x_jump=eta_jump*dvy_y[j];
double p_jump=-2*dvy_y[j]*eta_jump;
- double dzx_xx_jump=-dvx_yy[j]*eta_jump + p_jump;
+ double dp_x_jump=2*dvx_yy[j]*eta_jump;
+ double dzx_xx_jump=
+ -dvx_yy[j]*eta_jump + dp_x_jump;
dzx_xx_correction=-(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
if(x>middle)
dzx_xx_correction*=-1;
+ dzx_xx+=dzx_xx_correction;
// dzx_xx-=eta_jump*(zx_jump + dx*dzx_x_jump + dx*dx*dzx_xx_jump/2)/(h*h);
- // if(j==N/4)
- // std::cout << "dzx_xx "
- // << i << " "
- // << j << " "
- // << x << " "
- // << y << " "
- // // << zx[i+2][j]*exp(-log_etax[i+2][j]) << " "
- // // << zx[i+1][j]*exp(-log_etax[i+1][j]) << " "
- // // << zx[i][j]*exp(-log_etax[i][j]) << " "
- // // << zx[i-1][j]*exp(-log_etax[i-1][j]) << " "
- // // << zx[i-2][j]*exp(-log_etax[i-2][j]) << " "
- // // << zx[i-3][j]*exp(-log_etax[i-3][j]) << " "
- // // << 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]) << " "
- // // << p[i+1][j] << " "
- // // << p[i][j] << " "
- // // << p[i-1][j] << " "
- // // << p[i-2][j] << " "
+ if(j==Ny/4)
+ std::cout << "dzx_xx "
+ << i << " "
+ << j << " "
+ << x << " "
+ << dzx_xx << " "
+ // << dzx_xx_correction << " "
+ // << zx_jump/(h*h) << " "
+ // << dx*dzx_x_jump/(h*h) << " "
+ // << (dx*dx*dzx_xx_jump/2)/(h*h) << " "
- // // << 2*(zx[i+2][j]-zx[i+1][j])/h - p[i+1][j] << " "
- // // << 2*(zx[i+1][j]-zx[i][j])/h - p[i][j] << " "
- // // << 2*(zx[i][j]-zx[i-1][j])/h - p[i-1][j] << " "
- // // << 2*(zx[i-1][j]-zx[i-2][j])/h - p[i-2][j] << " "
+ // << -dvx_yy[j]*eta_jump << " "
+ // << dp_x_jump << " "
- // // << (zx[i+2][j] - zx[i+1][j])/h + (zy[i+1][j+1] - zy[i+1][j])/h << " "
- // // << (zx[i+1][j] - zx[i][j])/h + (zy[i][j+1] - zy[i][j])/h << " "
- // // << (zx[i][j] - zx[i-1][j])/h + (zy[i-1][j+1] - zy[i-1][j])/h << " "
- // // << (zx[i-1][j] - zx[i-2][j])/h + (zy[i-2][j+1] - zy[i-2][j])/h << " "
+ // << dzx_xx + zx_jump/(h*h) << " "
+ // << dzx_xx + zx_jump/(h*h) + dx*dzx_x_jump/(h*h) << " "
+ // << dzx_xx + zx_jump/(h*h) + dx*dzx_x_jump/(h*h) + (dx*dx*dzx_xx_jump/2)/(h*h) << " "
- // // << (zy[i+2][j+1]-zy[i+1][j+1])/h + (zx[i+2][j+1]-zx[i+2][j])/h << " "
- // // << (zy[i+1][j+1]-zy[i ][j+1])/h + (zx[i+1][j+1]-zx[i+1][j])/h << " "
- // // << (zy[i ][j+1]-zy[i-1][j+1])/h + (zx[i ][j+1]-zx[i ][j])/h << " "
- // // << (zy[i-1][j+1]-zy[i-2][j+1])/h + (zx[i-1][j+1]-zx[i-1][j])/h << " "
- // // << (zy[i-2][j+1]-zy[i-3][j+1])/h + (zx[i-2][j+1]-zx[i-2][j])/h << " "
+ // << dzx_xx + dzx_xx_correction << " "
-
- // // << (zy[i+2][j]-zy[i+1][j])/h << " "
- // // << (zx[i+2][j+1]-zx[i+2][j])/h << " "
- // // << (zy[i+1][j]-zy[i][j])/h << " "
- // // << (zx[i+1][j+1]-zx[i+1][j])/h << " "
- // // << (zy[i][j]-zy[i-1][j])/h << " "
- // // << (zx[i][j+1]-zx[i][j])/h << " "
- // // << (zy[i-1][j]-zy[i-2][j])/h << " "
- // // << (zx[i-1][j+1]-zx[i-1][j])/h << " "
- // // << (zy[i-2][j]-zy[i-3][j])/h << " "
- // // << (zx[i-2][j+1]-zx[i-2][j])/h << " "
-
- // // << zx[i+2][j] << " "
- // // << zx[i+1][j] << " "
- // // << zx[i][j] << " "
- // // << zx[i-1][j] << " "
- // // << zx[i-2][j] << " "
- // // << zy[i+2][j] << " "
- // // << zy[i+1][j] << " "
- // // << zy[i][j] << " "
- // // << zy[i-1][j] << " "
- // // << zy[i-2][j] << " "
-
- // // << dx << " "
- // // << dzx_xx << " "
- // // << dzx_xx_correction << " "
- // // << dzx_xx + dzx_xx_correction << " "
- // // << eta_jump*zx_jump/(h*h) << " "
- // // << dx*eta_jump*dzx_x_jump/(h*h) << " "
- // // << dx*dx*dzx_xx_jump/(2*h*h) << " "
- // // << (zx[i+2][j]-2*zx[i+1][j]+zx[i][j])/(h*h) << " "
- // // << (zx[i+3][j]-2*zx[i+2][j]+zx[i+1][j])/(h*h) << " "
- // // << (zx[i+4][j]-2*zx[i+3][j]+zx[i+2][j])/(h*h) << " "
- // // << (zx[i-3][j]-2*zx[i-2][j]+zx[i-1][j])/(h*h) << " "
- // << "\n";
+ << (zx[i+3][j]-2*zx[i+2][j]+zx[i+1][j])/(h*h) << " "
+ << (zx[i+2][j]-2*zx[i+1][j]+zx[i][j])/(h*h) << " "
+ << (zx[i+1][j]-2*zx[i][j]+zx[i-1][j])/(h*h) << " "
+ << (zx[i][j]-2*zx[i-1][j]+zx[i-2][j])/(h*h) << " "
+ << (zx[i-1][j]-2*zx[i-2][j]+zx[i-3][j])/(h*h) << " "
+ << "\n";
if(x+h/2>middle && x-h/2<middle)
{
dx=(x>middle) ? ((x-h/2)-middle)
: ((x+h/2)-middle);
- double dp_x_jump=2*dvx_yy[j]*eta_jump;
dp_x-=(p_jump + dx*dp_x_jump)/h;
-
- // if(j<3)
- // std::cout << "p "
- // << i << " "
- // << j << " "
- // << x << " "
- // << dp_x << " "
-
- // << (p[i+3][j]-p[i+2][j])/h << " "
- // << (p[i+2][j]-p[i+1][j])/h << " "
- // << (p[i+1][j]-p[i][j])/h << " "
- // << (p[i][j]-p[i-1][j])/h << " "
- // << (p[i-1][j]-p[i-2][j])/h << " "
- // << (p[i-2][j]-p[i-3][j])/h << " "
- // << "\n";
-
dzy_xy-=
eta_jump*(dvy_y[j] - dx*dvx_yy[j])/h;
- dzy_xy_correction=eta_jump*(dvy_y[j] - dx*dvx_yy[j])/h;
- if(j==N/4)
- std::cout << "dzy_xy "
- << i << " "
- << j << " "
- << x << " "
- << dzy_xy << " "
- // << dzy_xy_correction << " "
- // << dzy_xy - dzy_xy_correction << " "
-
- << ((zy[i+2][j+1] - zy[i+1][j+1])
- - (zy[i+2][j] - zy[i+1][j]))/(h*h) << " "
- << ((zy[i+1][j+1] - zy[i][j+1])
- - (zy[i+1][j] - zy[i][j]))/(h*h) << " "
- << ((zy[i][j+1] - zy[i-1][j+1])
- - (zy[i][j] - zy[i-1][j]))/(h*h) << " "
- << ((zy[i-1][j+1] - zy[i-2][j+1])
- - (zy[i-1][j] - zy[i-2][j]))/(h*h) << " "
- << "\n";
}
}
}
@@ -304,7 +215,7 @@ int main()
Resid_x[i][j]=Rx;
// zx[i][j]-=theta_mom*Rx/dRx_dzx;
- Resid_x[i][j]=dzy_xy;
+ // Resid_x[i][j]=dzx_xx;
}
}
}
@@ -314,18 +225,18 @@ int main()
{
std::stringstream ss;
ss << "zx_resid" << sweep;
- write_vtk(ss.str(),N+1,N,Resid_x);
+ write_vtk(ss.str(),Nx+1,Ny,Resid_x);
}
/* zy */
for(int rb=0;rb<2;++rb)
{
- for(int j=0;j<N+1;++j)
+ for(int j=0;j<Ny+1;++j)
{
int i_min=(j+rb)%2;
- for(int i=i_min;i<N;i+=2)
+ for(int i=i_min;i<Nx;i+=2)
{
- if(j!=0 && j!=N)
+ if(j!=0 && j!=Ny)
{
/* Do the finite difference parts */
@@ -345,7 +256,7 @@ int main()
dzy_xx=(-zy[i][j] + zy[i+1][j])/(h*h);
dzy_x=(zy[i+1][j] - zy[i][j])/(2*h);
}
- else if(i==N-1)
+ else if(i==Nx-1)
{
dzy_xx=(zy[i-1][j] - zy[i][j])/(h*h);
dzy_x=(zy[i][j] - zy[i-1][j])/(2*h);
@@ -430,13 +341,13 @@ int main()
{
std::stringstream ss;
ss << "zy_resid" << sweep;
- write_vtk(ss.str(),N,N,Resid_y);
+ write_vtk(ss.str(),Nx,Ny+1,Resid_y);
}
/* Pressure update */
- for(int j=0;j<N;++j)
- for(int i=0;i<N;++i)
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
{
double dzx_x=(zx[i+1][j] - zx[i][j])/h;
double dzy_y=(zy[i][j+1] - zy[i][j])/h;
@@ -517,14 +428,14 @@ int main()
{
std::stringstream ss;
ss << "p_resid" << sweep;
- write_vtk(ss.str(),N,N,Resid_p);
+ write_vtk(ss.str(),Nx,Ny,Resid_p);
}
/* Velocity Fix */
double max_x(0), max_y(0), max_p(0);
- for(int j=0;j<N;++j)
- for(int i=0;i<N;++i)
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
{
/* Fix vx */
if(i>0)
@@ -577,23 +488,23 @@ int main()
std::stringstream ss;
ss << "zx" << sweep;
- write_vtk(ss.str(),N+1,N,zx);
+ write_vtk(ss.str(),Nx+1,Ny,zx);
ss.str("");
ss << "zy" << sweep;
- write_vtk(ss.str(),N,N,zy);
+ write_vtk(ss.str(),Nx,Ny+1,zy);
ss.str("");
ss << "p" << sweep;
- write_vtk(ss.str(),N,N,p);
+ write_vtk(ss.str(),Nx,Ny,p);
- for(int i=0;i<N;++i)
- for(int j=0;j<N;++j)
+ for(int j=0;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
{
div[i][j]=(zx[i+1][j]/exp(log_etax[i+1][j])-zx[i][j]/exp(log_etax[i][j]))
+ (zy[i][j+1]/exp(log_etay[i][j+1])-zy[i][j]/exp(log_etay[i][j]));
}
ss.str("");
ss << "div" << sweep;
- write_vtk(ss.str(),N,N,div);
+ write_vtk(ss.str(),Nx,Ny,div);
}
if(fabs(max_x-max_x_previous)/max_x < tolerance
@@ -614,14 +525,14 @@ int main()
max_p_previous=max_p;
}
}
- for(int i=0;i<N+1;++i)
- for(int j=0;j<N+1;++j)
+ for(int i=0;i<Nx+1;++i)
+ for(int j=0;j<Ny+1;++j)
{
- if(j<N)
+ if(j<Ny)
zx[i][j]*=exp(-log_etax[i][j]);
- if(i<N)
+ if(i<Nx)
zy[i][j]*=exp(-log_etay[i][j]);
}
- write_vtk("vx",N+1,N,zx);
- write_vtk("vy",N,N,zy);
+ write_vtk("vx",Nx+1,Ny,zx);
+ write_vtk("vy",Nx,Ny+1,zy);
}
More information about the CIG-COMMITS
mailing list