[cig-commits] commit: Use Neumann boundaries on top and bottom for vx. Also start with zero for initial guess.
Mercurial
hg at geodynamics.org
Fri Mar 2 16:27:59 PST 2012
changeset: 73:fa26e2bf82f3
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 02 16:27:50 2012 -0800
files: compute_v_on_interface.cxx constants.hxx main.cxx
description:
Use Neumann boundaries on top and bottom for vx. Also start with zero for initial guess.
diff -r 91ce53ec8151 -r fa26e2bf82f3 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Fri Mar 02 12:19:55 2012 -0800
+++ b/compute_v_on_interface.cxx Fri Mar 02 16:27:50 2012 -0800
@@ -81,15 +81,15 @@ double compute_dv_yy(const double z[][ny
double dv_yy_p;
if(j==0)
{
- dv_yy_p=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
- - vel(z,log_eta,i+1,ny-1))/(h*h);
- // dv_yy_p=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
+ // dv_yy_p=(vel(z,log_eta,i+1,j+1) - 2*vel(z,log_eta,i+1,j)
+ // - vel(z,log_eta,i+1,ny-1))/(h*h);
+ dv_yy_p=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/(h*h);
}
else if(j==ny-1)
{
- dv_yy_p=(-vel(z,log_eta,i+1,0) - 2*vel(z,log_eta,i+1,j)
- + vel(z,log_eta,i+1,j-1))/(h*h);
- // dv_yy_p=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
+ // dv_yy_p=(-vel(z,log_eta,i+1,0) - 2*vel(z,log_eta,i+1,j)
+ // + vel(z,log_eta,i+1,j-1))/(h*h);
+ dv_yy_p=(-vel(z,log_eta,i+1,j) + vel(z,log_eta,i+1,j-1))/(h*h);
}
else
{
@@ -100,15 +100,15 @@ double compute_dv_yy(const double z[][ny
double dv_yy_m;
if(j==0)
{
- dv_yy_m=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
- - vel(z,log_eta,i,ny-1))/(h*h);
- // dv_yy_m=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
+ // dv_yy_m=(vel(z,log_eta,i,j+1) - 2*vel(z,log_eta,i,j)
+ // - vel(z,log_eta,i,ny-1))/(h*h);
+ dv_yy_m=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/(h*h);
}
else if(j==ny-1)
{
- dv_yy_m=(-vel(z,log_eta,i,0) - 2*vel(z,log_eta,i,j)
- + vel(z,log_eta,i,j-1))/(h*h);
- // dv_yy_m=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
+ // dv_yy_m=(-vel(z,log_eta,i,0) - 2*vel(z,log_eta,i,j)
+ // + vel(z,log_eta,i,j-1))/(h*h);
+ dv_yy_m=(-vel(z,log_eta,i,j) + vel(z,log_eta,i,j-1))/(h*h);
}
else
{
@@ -207,15 +207,6 @@ void compute_v_on_interface(const double
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();
@@ -226,15 +217,6 @@ void compute_v_on_interface(const double
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();
diff -r 91ce53ec8151 -r fa26e2bf82f3 constants.hxx
--- a/constants.hxx Fri Mar 02 12:19:55 2012 -0800
+++ b/constants.hxx Fri Mar 02 16:27:50 2012 -0800
@@ -1,8 +1,8 @@ const int N(256);
-const int N(256);
+const int N(128);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
-const double max_eta=1e6;
+const double max_eta=1e10;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
diff -r 91ce53ec8151 -r fa26e2bf82f3 main.cxx
--- a/main.cxx Fri Mar 02 12:19:55 2012 -0800
+++ b/main.cxx Fri Mar 02 16:27:50 2012 -0800
@@ -58,6 +58,18 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy);
+ for(int i=0;i<Nx+1;++i)
+ for(int j=0;j<Ny;++j)
+ zx[i][j]=0;
+
+ for(int i=0;i<Nx;++i)
+ for(int j=0;j<Ny+1;++j)
+ zy[i][j]=0;
+
+ for(int i=0;i<Nx;++i)
+ for(int j=0;j<Ny;++j)
+ p[i][j]=0;
+
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
for(int sweep=0;sweep<n_sweeps;++sweep)
@@ -88,19 +100,19 @@ int main()
double dzx_y, dzx_yy;
if(j==0)
{
- dzx_yy=
- (-zx[i][Ny-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
- dzx_y=(zx[i][j+1] + zx[i][Ny-1])/(2*h);
- // dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
- // dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
+ // dzx_yy=
+ // (-zx[i][Ny-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ // dzx_y=(zx[i][j+1] + zx[i][Ny-1])/(2*h);
+ 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==Ny-1)
{
- dzx_yy=
- (zx[i][j-1] - 2*zx[i][j] - zx[i][0])/(h*h);
- dzx_y=(-zx[i][0] - zx[i][j-1])/(2*h);
- // dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
- // dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
+ // dzx_yy=
+ // (zx[i][j-1] - 2*zx[i][j] - zx[i][0])/(h*h);
+ // dzx_y=(-zx[i][0] - zx[i][j-1])/(2*h);
+ dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
+ dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
}
else
{
More information about the CIG-COMMITS
mailing list