[cig-commits] commit: Seems to work for max_eta=10
Mercurial
hg at geodynamics.org
Wed Feb 22 13:30:54 PST 2012
changeset: 60:5ad7ab8e65a2
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 22 13:30:45 2012 -0800
files: compute_corrections.cxx constants.hxx main.cxx
description:
Seems to work for max_eta=10
diff -r 9d29fb7ccbb6 -r 5ad7ab8e65a2 compute_corrections.cxx
--- a/compute_corrections.cxx Wed Feb 22 05:23:50 2012 -0800
+++ b/compute_corrections.cxx Wed Feb 22 13:30:45 2012 -0800
@@ -20,26 +20,34 @@ void update_correction(const double &cur
old=current;
else
{
- if(current/old>1+theta)
- old*=1+theta;
- else if(current/old<1-theta)
- {
- if(current*old>1)
- old*=1-theta;
- else
- {
- double temp(old-(old-current)*theta);
- if(temp*old<0)
- {
- std::cout << "switch\n";
- old=0;
- }
- else
- old=temp;
- }
- }
- else
- old=current;
+ old+=theta*(current-old);
+
+ // double diff;
+ // if(current/old>1+theta)
+ // diff=theta*old;
+ // else if(current/old<1-theta)
+ // {
+ // if(current*old>1)
+ // diff=-theta*old;
+ // else
+ // {
+ // double temp(old-(old-current)*theta);
+ // if(temp*old<0)
+ // diff=-old;
+ // else
+ // diff=temp-old;
+ // }
+ // }
+ // else
+ // diff=current-old;
+
+ // double range(.1);
+ // if(diff>range)
+ // diff=range;
+ // else if(diff<-range)
+ // diff=-range;
+
+ // old+=diff;
}
}
diff -r 9d29fb7ccbb6 -r 5ad7ab8e65a2 constants.hxx
--- a/constants.hxx Wed Feb 22 05:23:50 2012 -0800
+++ b/constants.hxx Wed Feb 22 13:30:45 2012 -0800
@@ -1,15 +1,15 @@ const int N(512);
-const int N(512);
+const int N(128);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
-const double max_eta=2;
+const double max_eta=1e1;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
const double log_min_eta=std::log(min_eta);
// const double middle=(25 + 0.00000001)/64;
-// const double middle=(25 + 32/Nx + 0.00000001)/64;
-const double middle=(25 - 1/32.0)/64;
+const double middle=(25 + 0.00000001)/64 + 1/(2*Nx);
+// const double middle=(25 - 1/32.0)/64;
// const double middle=0.4;
const double h(1.0/N);
const double pi(atan(1.0)*4);
diff -r 9d29fb7ccbb6 -r 5ad7ab8e65a2 main.cxx
--- a/main.cxx Wed Feb 22 05:23:50 2012 -0800
+++ b/main.cxx Wed Feb 22 13:30:45 2012 -0800
@@ -38,8 +38,10 @@ int main()
int n_sweeps=10000000;
const double theta_mom=0.7;
const double theta_continuity=1.0;
- const double theta_correction=0.00001;
+ const double theta_correction=0.01;
const double tolerance=1.0e-6;
+
+ const int output_interval=100000000;
Model model(Model::solCx);
@@ -56,7 +58,8 @@ int main()
// if(iteration==max_iterations-1)
// n_sweeps=1;
double max_x_resid_previous(0), max_y_resid_previous(0), max_p_resid_previous(0);
- for(int sweep=0;sweep<n_sweeps;++sweep)
+ int sweep;
+ for(sweep=0;sweep<n_sweeps;++sweep)
{
// compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
@@ -179,14 +182,17 @@ int main()
}
}
- if(sweep%1000==0)
+ if(sweep%output_interval==0)
{
std::stringstream ss;
ss << "zx_resid" << sweep;
- write_vtk(ss.str(),Nx+1,N,Resid_x);
- ss.str("");
- ss << "Cx" << sweep;
- write_vtk(ss.str(),Nx+1,N,Cx);
+ // write_vtk(ss.str(),Nx+1,N,Resid_x);
+ if(sweep==0)
+ {
+ ss.str("");
+ ss << "Cx" << iteration;
+ write_vtk(ss.str(),Nx+1,N,Cx);
+ }
}
/* zy */
@@ -294,14 +300,17 @@ int main()
}
}
- if(sweep%1000==0)
+ if(sweep%output_interval==0)
{
std::stringstream ss;
ss << "zy_resid" << sweep;
- write_vtk(ss.str(),Nx,N,Resid_y);
- ss.str("");
- ss << "Cy" << sweep;
- write_vtk(ss.str(),Nx,N,Cy);
+ // write_vtk(ss.str(),Nx,N,Resid_y);
+ if(sweep==0)
+ {
+ ss.str("");
+ ss << "Cy" << iteration;
+ write_vtk(ss.str(),Nx,N,Cy);
+ }
}
/* Pressure update */
@@ -354,14 +363,17 @@ int main()
p[i][j]+=dp[i][j];
}
- if(sweep%1000==0)
+ if(sweep%output_interval==0)
{
std::stringstream ss;
ss << "p_resid" << sweep;
- write_vtk(ss.str(),Nx,N,Resid_p);
- ss.str("");
- ss << "Cp" << sweep;
- write_vtk(ss.str(),Nx,N,Cp);
+ // write_vtk(ss.str(),Nx,N,Resid_p);
+ if(sweep==0)
+ {
+ ss.str("");
+ ss << "Cp" << iteration;
+ write_vtk(ss.str(),Nx,N,Cp);
+ }
}
/* Velocity Fix */
@@ -430,27 +442,29 @@ int main()
<< max_y_resid << " "
<< max_p_resid << " "
<< "\n";
+ }
+ // if(sweep%output_interval==0)
+ // {
+ // std::stringstream ss;
+ // ss << "zx" << sweep;
+ // write_vtk(ss.str(),Nx+1,N,zx);
+ // ss.str("");
+ // ss << "zy" << sweep;
+ // write_vtk(ss.str(),Nx,N,zy);
+ // ss.str("");
+ // ss << "p" << sweep;
+ // write_vtk(ss.str(),Nx,N,p);
- std::stringstream ss;
- ss << "zx" << sweep;
- write_vtk(ss.str(),Nx+1,N,zx);
- ss.str("");
- ss << "zy" << sweep;
- write_vtk(ss.str(),Nx,N,zy);
- ss.str("");
- ss << "p" << sweep;
- write_vtk(ss.str(),Nx,N,p);
-
- 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(),Nx,N,div);
- }
+ // 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(),Nx,N,div);
+ // }
// if(fabs(max_x_resid-max_x_resid_previous)/max_x_resid < tolerance
// && fabs(max_y_resid-max_y_resid_previous)/max_y_resid < tolerance
@@ -462,9 +476,12 @@ int main()
std::cout << "Solved "
<< sweep << " "
<< tolerance << " "
- << fabs(max_x_resid-max_x_resid_previous)/max_x_resid << " "
- << fabs(max_y_resid-max_y_resid_previous)/max_y_resid << " "
- << fabs(max_p_resid-max_p_resid_previous)/max_p_resid << " "
+ << max_x_resid << " "
+ << max_y_resid << " "
+ << max_p_resid << " "
+ // << fabs(max_x_resid-max_x_resid_previous)/max_x_resid << " "
+ // << fabs(max_y_resid-max_y_resid_previous)/max_y_resid << " "
+ // << fabs(max_p_resid-max_p_resid_previous)/max_p_resid << " "
<< "\n";
break;
}
@@ -472,6 +489,8 @@ int main()
max_y_resid_previous=max_y_resid;
max_p_resid_previous=max_p_resid;
}
+ if(sweep==0)
+ break;
}
for(int i=0;i<Nx+1;++i)
for(int j=0;j<Ny+1;++j)
More information about the CIG-COMMITS
mailing list