[cig-commits] commit: Change to DIIM interpolation. Currently has analytic tests baked in for verification.
Mercurial
hg at geodynamics.org
Tue Feb 21 13:00:15 PST 2012
changeset: 54:f5cd216d64c9
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 21 11:21:25 2012 -0800
files: compute_v_on_interface.cxx constants.hxx main.cxx
description:
Change to DIIM interpolation. Currently has analytic tests baked in for verification.
diff -r 8f93ae19b917 -r f5cd216d64c9 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Tue Feb 21 11:18:56 2012 -0800
+++ b/compute_v_on_interface.cxx Tue Feb 21 11:21:25 2012 -0800
@@ -1,194 +1,208 @@
#include <algorithm>
#include <cmath>
+#include "constants.hxx"
+#include <limits>
#include <iostream>
-#include <iomanip>
-#include "constants.hxx"
-#include <map>
-#include <gsl/gsl_linalg.h>
-namespace {
- template<int Ny>
- void add_to_matrix(const int &n_terms, const int n, double z[][Ny],
- double log_eta[][Ny], const int i, const int j,
- const double &u0, const double &dx,
- const double &dy,
- double m[], double rhs[], const int sign=1)
- {
- if(j<0 || j>Ny-1)
- return;
- rhs[n]+=sign*z[i][j]*exp(-log_eta[i][j]);
- m[0+n*n_terms]+=sign*u0;
- m[1+n*n_terms]+=sign*dy;
- m[2+n*n_terms]+=sign*dy*dy/2;
- m[3+n*n_terms]+=sign*dx;
- m[4+n*n_terms]+=sign*dx*dy;
- m[5+n*n_terms]+=sign*dx*dx/2;
- }
-}
+double analytic(const double x, const double y, const bool return_ux)
+{
+ double ux(1.1), ux_xp(1.2), ux_xxp(1.3), ux_xyp(1.4),
+ ux_xm(1.5), ux_xxm(1.6), ux_xym(1.7),
+ uy(2.1), uy_xp(2.2), uy_xxp(2.3), uy_xyp(2.4),
+ uy_xm(2.5), uy_xxm(2.6), uy_xym(2.7);
+ double ux_y(-(max_eta*uy_xp - min_eta*uy_xm)/eta_jump),
+ ux_yy(-(max_eta*uy_xyp - min_eta*uy_xym)/eta_jump),
+ uy_y(-(max_eta*ux_xp - min_eta*ux_xm)/eta_jump),
+ uy_yy(-(max_eta*ux_xyp - min_eta*ux_xym)/eta_jump);
-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;
-}
+ // std::cout << "u y "
+ // << ux << " "
+ // << ux_y << " "
+ // << ux_yy << " "
+ // << uy << " "
+ // << uy_y << " "
+ // << uy_yy << " "
+ // << "\n";
+ // exit(0);
-template<int Ny>
-void interpolate_to_surface(double z[][Ny], double log_eta[][Ny],
- const int &i, const int &j,
- const double &delta_x, const bool &aligned,
- double &v0, double &dv_y, double &dv_yy)
-{
- const int n_terms(6);
- /* Format of m is
-
- u0, u,y u,yy um,x um,xy um,xx up,x up,xy up,xx
-
- where y is parallel to the surface and x is
- perpendicular */
-
- double m[n_terms*n_terms];
- double rhs[n_terms];
-
- std::fill(m,m+n_terms*n_terms,0.0);
- std::fill(rhs,rhs+n_terms,0.0);
-
- double dx(delta_x);
- int i_center(i), di(1);
- if(dx>h/2)
+ if(return_ux)
{
- di=-1;
- i_center++;
- dx=h-dx;
- }
-
- if(aligned)
- {
- /* Ordering is
- i_center,j
- i_center,j+1
- i_center,j-1
- i_center-1,j
- i_center-2,j
- i_center-1,j+1 - i_center+1,j-1
-
- i_center+1,j
- i_center+2,j
- i_center+1,j+1 - i_center-1,j-1
- */
-
- add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,di*dx,0,m,rhs);
- add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,di*dx,h,m,rhs);
- add_to_matrix(n_terms,2,z,log_eta,i_center,j-1,1,di*dx,-h,m,rhs);
- add_to_matrix(n_terms,3,z,log_eta,i_center-di,j,1,di*(dx+h),0,m,rhs);
- add_to_matrix(n_terms,4,z,log_eta,i_center-2*di,j,1,di*(dx+2*h),
- 0,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-di,j+1,1,di*(dx+h),h,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-di,j-1,1,di*(dx+h),
- -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);
-
- // }
-
- /* At the j==0 boundary, all of the j-1 terms are not applied.
- Instead, for n=2, apply dvx/dy=0. */
- if(j==0)
- {
- m[1+2*n_terms]=1;
- m[2+2*n_terms]=-h/2;
- }
- else if(j==Ny-1)
- {
- m[1+1*n_terms]=1;
- m[2+1*n_terms]=h/2;
- }
+ if(x>0)
+ return ux + x*ux_xp + x*x*ux_xxp/2 + x*y*ux_xyp
+ + y*ux_y + y*y*ux_yy/2;
+ else
+ return ux + x*ux_xm + x*x*ux_xxm/2 + x*y*ux_xym
+ + y*ux_y + y*y*ux_yy/2;
}
else
{
- /* Ordering is
- i_center,j
- i_center,j+1
- i_center-1,j
- i_center-1,j+1
- i_center,j+2 + i_center,j-1
- i_center-2,j + i_center-2,j+1
-
- i_center+1,j
- i_center+1,j+1
- i_center+2,j + i_center+2,j+1
- */
-
- add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,-di*dx,-h/2,m,rhs);
- add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,-di*dx,h/2,m,rhs);
- add_to_matrix(n_terms,2,z,log_eta,i_center-di,j,1,-di*(dx+h),
- -h/2,m,rhs);
- add_to_matrix(n_terms,3,z,log_eta,i_center-di,j+1,1,-di*(dx+h),
- h/2,m,rhs);
- add_to_matrix(n_terms,4,z,log_eta,i_center,j+2,1,-di*dx,3*h/2,m,rhs);
- add_to_matrix(n_terms,4,z,log_eta,i_center,j-1,1,-di*dx,-3*h/2,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j,1,-di*(dx+2*h),
- -h/2,m,rhs);
- add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j+1,1,-di*(dx+2*h),
- h/2,m,rhs);
+ if(x>0)
+ return uy + x*uy_xp + x*x*uy_xxp/2 + x*y*uy_xyp
+ + y*uy_y + y*y*uy_yy/2;
+ else
+ return uy + x*uy_xm + x*x*uy_xxm/2 + x*y*uy_xym
+ + y*uy_y + y*y*uy_yy/2;
}
- // if(j>0)
- // for(int a=0;a<n_terms;++a)
- // {
- // std::cout << a << " "
- // << std::setw(15)
- // << std::setprecision(10)
- // << rhs[a] << " ";
- // for(int b=0;b<n_terms;++b)
- // std::cout << std::setw(15)
- // << std::setprecision(10)
- // << m[b+n_terms*a] << " ";
- // std::cout << "\n";
- // }
-
- gsl_matrix_view mv=gsl_matrix_view_array(m,n_terms,n_terms);
- gsl_vector_view rhsv=gsl_vector_view_array(rhs,n_terms);
- gsl_vector *result=gsl_vector_alloc(n_terms);
- int s;
- 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<n_terms;++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);
- dv_yy=gsl_vector_get(result,2);
- gsl_vector_free(result);
- gsl_permutation_free(perm);
}
+
+template<int ny>
+double vel(double z[][ny],double log_eta[][ny], const int i, const int j)
+{
+ if(j>ny/8-4 && j<ny/8+4)
+ {
+ double x,y;
+ if(ny%2==0)
+ {
+ x=i*h-middle;
+ y=(j+0.5)*h-1.0/8;
+ // y=j*h-1.0/8;
+ return analytic(x,y,true);
+ }
+ else
+ {
+ x=(i+0.5)*h-middle;
+ y=j*h-1.0/8;
+ // y=(j-0.5)*h-1.0/8;
+ return analytic(x,y,false);
+ }
+ }
+
+ return z[i][j]*std::exp(-log_eta[i][j]);
+}
+
+template<int ny>
+double compute_dv_yy(double z[][ny],
+ double log_eta[][ny], const double &dx,
+ const int i, const int j)
+{
+ double dv_yy_p1=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
+ + vel(z,log_eta,i+1,j-1))/(h*h);
+ double dv_yy_p2=(vel(z,log_eta,i+2,j+1) - 2*vel(z,log_eta,i+2,j)
+ + vel(z,log_eta,i+2,j-1))/(h*h);
+ if(j==0)
+ {
+ dv_yy_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
+ dv_yy_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/(h*h);
+ }
+ else if(j==ny-1)
+ {
+ dv_yy_p1=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
+ dv_yy_p2=(-vel(z,log_eta,i+2,j) + vel(z,log_eta,i+2,j-1))/(h*h);
+ }
+ const double dv_yy_p=(1-dx)*dv_yy_p1 + dx*dv_yy_p2;
+
+ double dv_yy_m0=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
+ + vel(z,log_eta,i,j-1))/(h*h);
+ double dv_yy_m1=(vel(z,log_eta,i-1,j+1) - 2*vel(z,log_eta,i-1,j)
+ + vel(z,log_eta,i-1,j-1))/(h*h);
+ if(j==0)
+ {
+ dv_yy_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
+ dv_yy_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/(h*h);
+ }
+ else if(j==ny-1)
+ {
+ dv_yy_m0=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
+ dv_yy_m1=(-vel(z,log_eta,i-1,j) + vel(z,log_eta,i-1,j-1))/(h*h);
+ }
+ double dv_yy_m=(1-dx)*dv_yy_m1 + dx*dv_yy_m0;
+
+ // if(j==ny/8)
+ // std::cout << "dvy_yy "
+ // << ny << " "
+ // << dv_yy_m << " "
+ // << dv_yy_m0 << " "
+ // << dv_yy_m1 << " "
+ // << dv_yy_p << " "
+ // << dv_yy_p1 << " "
+ // << dv_yy_p2 << " "
+ // << (dv_yy_p + dv_yy_m)/2 << " "
+ // << "\n";
+
+ return (dv_yy_p + dv_yy_m)/2;
+}
+
+template<int ny>
+double compute_dv_y(double z[][ny],
+ double log_eta[][ny],
+ const double &jump,
+ const double &dx,
+ const int i, const int j)
+{
+ double dv_y_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/h;
+ double dv_y_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/h;
+ double dv_y_p=(1-dx)*dv_y_p1 + dx*dv_y_p2;
+
+ double dv_y_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/h;
+ double dv_y_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/h;
+ double dv_y_m=(1-dx)*dv_y_m1 + dx*dv_y_m0;
+
+ if(j==ny/8 && ny%2==1)
+ std::cout << "dvy_y "
+ << ny << " "
+ << i*h-middle << " "
+ << (j+0.5)*h-1.0/8 << " "
+ << (j+1.5)*h-1.0/8 << " "
+ << vel(z,log_eta,i+1,j+1) << " "
+ << vel(z,log_eta,i+1,j) << " "
+ << vel(z,log_eta,i+2,j+1) << " "
+ << vel(z,log_eta,i+2,j) << " "
+ << dv_y_m << " "
+ << dv_y_m0 << " "
+ << dv_y_m1 << " "
+ << dv_y_p << " "
+ << dv_y_p1 << " "
+ << dv_y_p2 << " "
+ << (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
+ /(max_eta + min_eta) << " "
+ << "\n";
+
+ return (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
+ /(max_eta + min_eta);
+}
+
+template<int ny>
+double compute_v(double z[][ny],
+ double log_eta[][ny],
+ const double &jump,
+ const double &dx,
+ const int i, const int j)
+{
+ double v_p1(vel(z,log_eta,i+1,j)), v_p2(vel(z,log_eta,i+2,j)),
+ v_p3(vel(z,log_eta,i+3,j));
+ double v_p((1-dx)*v_p1 + dx*v_p2), v_pp((1-dx)*v_p2 + dx*v_p3);
+
+ double v_m0(vel(z,log_eta,i,j)), v_m1(vel(z,log_eta,i-1,j)),
+ v_m2(vel(z,log_eta,i-2,j));
+ double v_m((1-dx)*v_m1 + dx*v_m0), v_mm((1-dx)*v_m2 + dx*v_m1);
+
+ return (4*(max_eta*v_p + min_eta*v_m) - (max_eta*v_pp + min_eta*v_mm)
+ - 2*h*jump)
+ /(3*(max_eta+min_eta));
+}
+
+/* For now, hard code the solCx answer */
+double dzy_yx_jump(const double &dvx_yy)
+{
+ return -eta_jump*dvx_yy;
+}
+
+double dzx_x_jump(const double &dvy_y)
+{
+ return -eta_jump*dvy_y;
+}
+
+double dzx_yx_jump(const double &dvy_yy)
+{
+ return -eta_jump*dvy_yy;
+}
+
+double dzy_x_jump(const double &dvx_y)
+{
+ return -eta_jump*dvx_y;
+}
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],
@@ -198,35 +212,44 @@ void compute_v_on_interface(double zx[Nx
double &dvx_y, double &dvy_y,
double &dvx_yy, double &dvy_yy)
{
- int jx,jy;
+ int ix(x/h), iy(x/h - 0.5);
+ double dx((x-ix*h)/h), dy((x-iy*h)/h-0.5);
if(xyz==0)
{
- jx=jy=j;
+ dvx_yy=compute_dv_yy(zx,log_etax,dx,ix,j);
+ dvy_y=compute_dv_y(zy,log_etay,dzy_yx_jump(dvx_yy),dy,iy,j);
+ vx=compute_v(zx,log_etax,dzx_x_jump(dvy_y),dx,ix,j);
+
+ if(j==Ny/8)
+ std::cout << "compute x "
+ << x << " "
+ << j << " "
+ << vx << " "
+ << dvy_y << " "
+ << dvx_yy << " "
+ << "\n";
+
+ vy=std::numeric_limits<double>::infinity();
+ dvx_y=std::numeric_limits<double>::infinity();
+ dvy_yy=std::numeric_limits<double>::infinity();
}
else
{
- jx=j-1;
- jy=j;
+ dvy_yy=compute_dv_yy(zy,log_etay,dy,iy,j);
+ dvx_y=compute_dv_y(zx,log_etax,dzx_yx_jump(dvy_yy),dx,ix,j-1);
+ vy=compute_v(zy,log_etay,dzy_x_jump(dvx_y),dy,iy,j);
+
+ if(j==Ny/8)
+ std::cout << "compute y "
+ << x << " "
+ << j << " "
+ << vy << " "
+ << dvx_y << " "
+ << dvy_yy << " "
+ << "\n";
+
+ vx=std::numeric_limits<double>::infinity();
+ dvy_y=std::numeric_limits<double>::infinity();
+ dvx_yy=std::numeric_limits<double>::infinity();
}
-
- // std::cout << "compute v "
- // << x << " "
- // << j << " "
- // << xyz << " "
- // << "\n";
-
- /* vx */
- int ix(x/h);
- double dx_x(x-ix*h);
- interpolate_to_surface(zx,log_etax,ix,jx,dx_x,xyz==0,vx,dvx_y,dvx_yy);
-
- // std::cout
- // << vx << " "
- // << "\n";
-
-
- /* vy */
- int iy=x/h-0.5;
- double dx_y=x-((iy+0.5)*h);
- interpolate_to_surface(zy,log_etay,iy,jy,dx_y,xyz!=0,vy,dvy_y,dvy_yy);
}
diff -r 8f93ae19b917 -r f5cd216d64c9 constants.hxx
--- a/constants.hxx Tue Feb 21 11:18:56 2012 -0800
+++ b/constants.hxx Tue Feb 21 11:21:25 2012 -0800
@@ -1,8 +1,8 @@ const int NN(128);
-const int NN(128);
-const int Nx(NN);
-const int Ny(NN);
+const int N(256);
+const int Nx(N);
+const int Ny(N);
const double min_eta=1;
-const double max_eta=1e3;
+const double max_eta=1e4;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
@@ -11,7 +11,7 @@ const double log_min_eta=std::log(min_et
// const double middle=(25 + 32/Nx + 0.00000001)/64;
// const double middle=(25 - 32/2048.0)/64;
const double middle=0.4;
-const double h(1.0/Nx);
+const double h(1.0/N);
const double pi(atan(1.0)*4);
const double r2_inclusion=0.1;
diff -r 8f93ae19b917 -r f5cd216d64c9 main.cxx
--- a/main.cxx Tue Feb 21 11:18:56 2012 -0800
+++ b/main.cxx Tue Feb 21 11:21:25 2012 -0800
@@ -32,7 +32,7 @@ int main()
/* Initial conditions */
const int max_iterations=1;
- int n_sweeps=10000000;
+ int n_sweeps=1;
const double theta_mom=0.7/eta_jump;
const double theta_c=1.0/eta_jump;
const double tolerance=1.0e-6;
@@ -171,17 +171,37 @@ int main()
// << middle << " "
// // << zx_jump << " "
// // << dzx_x_jump << " "
- // << dzx_xx_jump << " "
- // << pi*pi*zx[i+1][j]*exp(-log_etax[i+1][j])*eta_jump << " "
- // << pi*pi*zx[i][j]*exp(-log_etax[i][j])*eta_jump << " "
- // << pi*pi*zx[i-1][j]*exp(-log_etax[i-1][j])*eta_jump << " "
+ // // << dzx_xx_jump << " "
+ // // << pi*pi*zx[i+1][j]*exp(-log_etax[i+1][j])*eta_jump << " "
+ // // << pi*pi*zx[i][j]*exp(-log_etax[i][j])*eta_jump << " "
+ // // << pi*pi*zx[i-1][j]*exp(-log_etax[i-1][j])*eta_jump << " "
// // << dzx_xx << " "
- // << (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) << " "
- // << dzx_xx+dzx_xx_correction << " "
+ // // << (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) << " "
+ // // << dzx_xx+dzx_xx_correction << " "
+
+ // // << vx << " "
+ // // << exp(-log_etax[i+2][j])*zx[i+2][j] << " "
+ // // << exp(-log_etax[i+1][j])*zx[i+1][j] << " "
+ // // << exp(-log_etax[i][j])*zx[i][j] << " "
+ // // << exp(-log_etax[i-1][j])*zx[i-1][j] << " "
+ // // << exp(-log_etax[i-2][j])*zx[i-2][j] << " "
+
+ // // << dvy_y << " "
+ // // << exp(-log_etay[i+2][j])*(zy[i+2][j+1]-zy[i+2][j])/(h) << " "
+ // // << exp(-log_etay[i+1][j])*(zy[i+1][j+1]-zy[i+1][j])/(h) << " "
+ // // << exp(-log_etay[i][j])*(zy[i][j+1]-zy[i][j])/(h) << " "
+ // // << exp(-log_etay[i-1][j])*(zy[i-1][j+1]-zy[i-1][j])/(h) << " "
+
+
+ // // << dvx_yy << " "
+ // // << exp(-log_etax[i+2][j])*(zx[i+2][j+1]-2*zx[i+2][j]+zx[i+2][j-1])/(h*h) << " "
+ // // << exp(-log_etax[i+1][j])*(zx[i+1][j+1]-2*zx[i+1][j]+zx[i+1][j-1])/(h*h) << " "
+ // // << exp(-log_etax[i][j])*(zx[i][j+1]-2*zx[i][j]+zx[i][j-1])/(h*h) << " "
+ // // << exp(-log_etax[i-1][j])*(zx[i-1][j+1]-2*zx[i-1][j]+zx[i-1][j-1])/(h*h) << " "
// << "\n";
@@ -191,6 +211,20 @@ int main()
{
dx=(x>middle) ? ((x-h/2)-middle)
: ((x+h/2)-middle);
+
+ // if(j==Ny/8)
+ // std::cout << "dp "
+ // << i << " "
+ // << j << " "
+ // << x << " "
+ // << middle << " "
+ // << dp_x-(p_jump + dx*dp_x_jump)/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";
dp_x-=(p_jump + dx*dp_x_jump)/h;
@@ -210,7 +244,7 @@ int main()
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
Resid_x[i][j]=Rx;
- zx[i][j]-=theta*Rx/dRx_dzx;
+ // zx[i][j]-=theta*Rx/dRx_dzx;
}
}
}
@@ -220,7 +254,7 @@ int main()
{
std::stringstream ss;
ss << "zx_resid" << sweep;
- write_vtk(ss.str(),Nx+1,NN,Resid_x);
+ write_vtk(ss.str(),Nx+1,N,Resid_x);
}
/* zy */
@@ -332,33 +366,55 @@ int main()
: ((x+h/2)-middle);
double dzx_yx_correction=
eta_jump*(dvx_y - dx*dvy_yy)/h;
+ if(j==Ny/8)
+ std::cout << "dzx_yx "
+ << i << " "
+ << j << " "
+ << dx << " "
+ << middle << " "
+ << x << " "
+ // << ((x-h/2)-middle)/dx << " "
+ // << ((x+h/2)-middle)/dx << " "
+ // << dzx_yx << " "
+ // << dzx_yx_correction << " "
+ // << eta_jump*dvx_y[j]/h << " "
+ // << eta_jump*dx*dvy_yy[j]/h << " "
+ // << dzx_yx << " "
+ // << eta_jump*(dvx_y)/h << " "
+ // << eta_jump*(-dx*dvy_yy)/h << " "
+ << dzx_yx - dzx_yx_correction << " "
+
+ << ((zx[i+3][j] - zx[i+3][j-1])
+ - (zx[i+2][j] - zx[i+2][j-1]))/(h*h) << " "
+ << ((zx[i+2][j] - zx[i+2][j-1])
+ - (zx[i+1][j] - zx[i+1][j-1]))/(h*h) << " "
+ << ((zx[i+1][j] - zx[i+1][j-1])
+ - (zx[i][j] - zx[i][j-1]))/(h*h) << " "
+ << ((zx[i][j] - zx[i][j-1])
+ - (zx[i-1][j] - zx[i-1][j-1]))/(h*h) << " "
+ << ((zx[i-1][j] - zx[i-1][j-1])
+ - (zx[i-2][j] - zx[i-2][j-1]))/(h*h) << " "
+ << ((zx[i-2][j] - zx[i-2][j-1])
+ - (zx[i-3][j] - zx[i-3][j-1]))/(h*h) << " "
+
+ << dvx_y << " "
+ << exp(-log_etax[i+2][j])*(zx[i+2][j]-zx[i+2][j-1])/h << " "
+ << exp(-log_etax[i+1][j])*(zx[i+1][j]-zx[i+1][j-1])/h << " "
+ << exp(-log_etax[i][j])*(zx[i][j]-zx[i][j-1])/h << " "
+ << exp(-log_etax[i-1][j])*(zx[i-1][j]-zx[i-1][j-1])/h << " "
+ << exp(-log_etax[i-2][j])*(zx[i-2][j]-zx[i-2][j-1])/h << " "
+
+
+ // << dvy_yy << " "
+ // << exp(-log_etay[i+2][j])*(zy[i+2][j+1]-2*zy[i+2][j]+zy[i+2][j-1])/(h*h) << " "
+ // << exp(-log_etay[i+1][j])*(zy[i+1][j+1]-2*zy[i+1][j]+zy[i+1][j-1])/(h*h) << " "
+ // << exp(-log_etay[i][j])*(zy[i][j+1]-2*zy[i][j]+zy[i][j-1])/(h*h) << " "
+ // << exp(-log_etay[i-1][j])*(zy[i-1][j+1]-2*zy[i-1][j]+zy[i-1][j-1])/(h*h) << " "
+ // << exp(-log_etay[i-2][j])*(zy[i-2][j+1]-2*zy[i-2][j]+zy[i-2][j-1])/(h*h) << " "
+
+
+ << "\n";
dzx_yx-=dzx_yx_correction;
- // if(j==Ny/8)
- // std::cout << "dzx_yx "
- // << i << " "
- // << j << " "
- // << middle << " "
- // << x << " "
- // // << ((x-h/2)-middle)/dx << " "
- // // << ((x+h/2)-middle)/dx << " "
- // // << dzx_yx << " "
- // // << dzx_yx_correction << " "
- // // << eta_jump*dvx_y[j]/h << " "
- // // << eta_jump*dx*dvy_yy[j]/h << " "
- // << dzx_yx << " "
- // << ((zx[i+3][j] - zx[i+3][j-1])
- // - (zx[i+2][j] - zx[i+2][j-1]))/(h*h) << " "
- // << ((zx[i+2][j] - zx[i+2][j-1])
- // - (zx[i+1][j] - zx[i+1][j-1]))/(h*h) << " "
- // << ((zx[i+1][j] - zx[i+1][j-1])
- // - (zx[i][j] - zx[i][j-1]))/(h*h) << " "
- // << ((zx[i][j] - zx[i][j-1])
- // - (zx[i-1][j] - zx[i-1][j-1]))/(h*h) << " "
- // << ((zx[i-1][j] - zx[i-1][j-1])
- // - (zx[i-2][j] - zx[i-2][j-1]))/(h*h) << " "
- // << ((zx[i-2][j] - zx[i-2][j-1])
- // - (zx[i-3][j] - zx[i-3][j-1]))/(h*h) << " "
- // << "\n";
}
}
@@ -393,7 +449,7 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
Resid_y[i][j]=Ry;
- zy[i][j]-=theta*Ry/dRy_dzy;
+ // zy[i][j]-=theta*Ry/dRy_dzy;
}
}
}
@@ -403,7 +459,7 @@ int main()
{
std::stringstream ss;
ss << "zy_resid" << sweep;
- write_vtk(ss.str(),Nx,NN,Resid_y);
+ write_vtk(ss.str(),Nx,N,Resid_y);
}
/* Pressure update */
@@ -477,14 +533,14 @@ int main()
Resid_p[i][j]=Rc;
dp[i][j]=-theta*Rc/dRc_dp;
- p[i][j]+=dp[i][j];
+ // p[i][j]+=dp[i][j];
}
if(sweep%1000==0 || sweep>236000)
{
std::stringstream ss;
ss << "p_resid" << sweep;
- write_vtk(ss.str(),Nx,NN,Resid_p);
+ write_vtk(ss.str(),Nx,N,Resid_p);
}
/* Velocity Fix */
@@ -514,7 +570,7 @@ int main()
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
- zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+ // zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
max_x_resid=std::max(std::fabs(Resid_x[i][j]),max_x_resid);
}
@@ -538,7 +594,7 @@ int main()
double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
- zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+ // zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
max_y_resid=std::max(fabs(Resid_y[i][j]),max_y_resid);
}
@@ -556,13 +612,13 @@ int main()
std::stringstream ss;
ss << "zx" << sweep;
- write_vtk(ss.str(),Nx+1,NN,zx);
+ write_vtk(ss.str(),Nx+1,N,zx);
ss.str("");
ss << "zy" << sweep;
- write_vtk(ss.str(),Nx,NN,zy);
+ write_vtk(ss.str(),Nx,N,zy);
ss.str("");
ss << "p" << sweep;
- write_vtk(ss.str(),Nx,NN,p);
+ write_vtk(ss.str(),Nx,N,p);
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i)
@@ -572,7 +628,7 @@ int main()
}
ss.str("");
ss << "div" << sweep;
- write_vtk(ss.str(),Nx,NN,div);
+ write_vtk(ss.str(),Nx,N,div);
}
if(fabs(max_x_resid-max_x_resid_previous)/max_x_resid < tolerance
@@ -601,6 +657,6 @@ int main()
if(i<Nx)
zy[i][j]*=exp(-log_etay[i][j]);
}
- write_vtk("vx",Nx+1,NN,zx);
- write_vtk("vy",Nx,NN,zy);
+ write_vtk("vx",Nx+1,N,zx);
+ write_vtk("vy",Nx,N,zy);
}
More information about the CIG-COMMITS
mailing list