[cig-commits] commit: dzy_xx correction works. Also fixed some issues with write_vtk offsets and minor cleanups
Mercurial
hg at geodynamics.org
Fri Feb 10 16:00:42 PST 2012
changeset: 41:37df131d70b5
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 08 12:20:45 2012 -0800
files: main.cxx
description:
dzy_xx correction works. Also fixed some issues with write_vtk offsets and minor cleanups
diff -r fa1aec3271c8 -r 37df131d70b5 main.cxx
--- a/main.cxx Wed Feb 08 12:19:45 2012 -0800
+++ b/main.cxx Wed Feb 08 12:20:45 2012 -0800
@@ -42,9 +42,6 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
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;
for(int iteration=0;iteration<max_iterations;++iteration)
{
@@ -156,7 +153,7 @@ int main()
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=
+ double dzx_xx_correction=
-(zx_jump + dx*dzx_x_jump
+ dx*dx*dzx_xx_jump/2)/(h*h);
if(x>middle)
@@ -196,7 +193,7 @@ int main()
{
std::stringstream ss;
ss << "zx_resid" << sweep;
- write_vtk(ss.str(),Nx+1,Ny,Resid_x);
+ write_vtk(ss.str(),Nx+1,NN,Resid_x);
}
/* zy */
@@ -282,16 +279,71 @@ int main()
: (middle-(x+h));
double zy_jump=eta_jump*vy[j];
double dzy_x_jump=-eta_jump*dvx_y[j];
- double dzy_xx_jump=-eta_jump*dvy_yy[j];
- dzy_xx_correction=(zy_jump + dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
- // dzy_xx-=(zy_jump + dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
+ double dzy_xx_jump=-3*eta_jump*dvy_yy[j];
+ double dzy_xx_correction=
+ -(zy_jump - dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
+ if(x>middle)
+ dzy_xx_correction*=-1;
+ dzy_xx+=dzy_xx_correction;
+
+ if(j==Ny/8)
+ std::cout << "dzy_xx "
+ << i << " "
+ << j << " "
+ << x << " "
+ // << middle << " "
+ // << h << " "
+ << dx << " "
+ << dzy_xx << " "
+ // << dzy_xx + dzy_xx_correction << " "
+ // << dzy_xx + zy_jump/(h*h) - dx*dzy_x_jump/(h*h)
+ // + dx*dx*dzy_xx_jump/(2*h*h) << " "
+ // << zy_jump/(h*h) << " "
+ // << dx*dzy_x_jump/(h*h) << " "
+ << dx*dx*dzy_xx_jump/(2*h*h) << " "
+ // << dzy_xx + dzy_xx_correction << " "
+ << (zy[i+2][j] - 2*zy[i+3][j] + zy[i+4][j])/(h*h) << " "
+ << (zy[i+1][j] - 2*zy[i+2][j] + zy[i+3][j])/(h*h) << " "
+ << (zy[i][j] - 2*zy[i+1][j] + zy[i+2][j])/(h*h) << " "
+ << (zy[i-1][j] - 2*zy[i][j] + zy[i+1][j])/(h*h) << " "
+ << (zy[i-2][j] - 2*zy[i-1][j] + zy[i][j])/(h*h) << " "
+ << (zy[i-3][j] - 2*zy[i-2][j] + zy[i-1][j])/(h*h) << " "
+ << (zy[i-4][j] - 2*zy[i-3][j] + zy[i-2][j])/(h*h) << " "
+ << "\n";
+
if(x+h/2>middle && x-h/2<middle)
{
dx=(x>middle) ? ((x-h/2)-middle)
: ((x+h/2)-middle);
- dzx_yx_correction=eta_jump*(dvx_y[j] + dx*dvy_yy[j])/h;
+ double dzx_yx_correction=eta_jump*(dvx_y[j] + dx*dvy_yy[j])/h;
// dzx_yx-=eta_jump*(dvx_y[j] + dx*dvy_yy[j])/h;
+ // if(j==Ny/4)
+ // 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 - 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) << " "
+ // << "\n";
}
}
}
@@ -308,6 +360,8 @@ int main()
Resid_y[i][j]=Ry;
// zy[i][j]-=theta_mom*Ry/dRy_dzy;
+
+ Resid_y[i][j]=dzy_xx;
}
}
}
@@ -317,7 +371,7 @@ int main()
{
std::stringstream ss;
ss << "zy_resid" << sweep;
- write_vtk(ss.str(),Nx,Ny+1,Resid_y);
+ write_vtk(ss.str(),Nx,NN,Resid_y);
}
/* Pressure update */
@@ -353,7 +407,7 @@ int main()
double zx_jump=vx[j+1];
double dzx_x_jump=dvy_y[j];
- dzx_x_correction=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
+ double dzx_x_correction=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
// dzx_x-=eta_jump*(zx_jump + dx*dzx_x_jump)/h;
// if(j==N/4)
@@ -404,7 +458,7 @@ int main()
{
std::stringstream ss;
ss << "p_resid" << sweep;
- write_vtk(ss.str(),Nx,Ny,Resid_p);
+ write_vtk(ss.str(),Nx,NN,Resid_p);
}
/* Velocity Fix */
@@ -464,13 +518,13 @@ int main()
std::stringstream ss;
ss << "zx" << sweep;
- write_vtk(ss.str(),Nx+1,Ny,zx);
+ write_vtk(ss.str(),Nx+1,NN,zx);
ss.str("");
ss << "zy" << sweep;
- write_vtk(ss.str(),Nx,Ny+1,zy);
+ write_vtk(ss.str(),Nx,NN,zy);
ss.str("");
ss << "p" << sweep;
- write_vtk(ss.str(),Nx,Ny,p);
+ write_vtk(ss.str(),Nx,NN,p);
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i)
@@ -480,7 +534,7 @@ int main()
}
ss.str("");
ss << "div" << sweep;
- write_vtk(ss.str(),Nx,Ny,div);
+ write_vtk(ss.str(),Nx,NN,div);
}
if(fabs(max_x-max_x_previous)/max_x < tolerance
@@ -509,6 +563,6 @@ int main()
if(i<Nx)
zy[i][j]*=exp(-log_etay[i][j]);
}
- write_vtk("vx",Nx+1,Ny,zx);
- write_vtk("vy",Nx,Ny+1,zy);
+ write_vtk("vx",Nx+1,NN,zx);
+ write_vtk("vy",Nx,NN,zy);
}
More information about the CIG-COMMITS
mailing list