[cig-commits] commit: Compute dC_dv and test it
Mercurial
hg at geodynamics.org
Wed Feb 29 14:54:22 PST 2012
changeset: 61:674bdfcc3cb8
user: Walter Landry <wlandry at caltech.edu>
date: Tue Feb 28 13:58:32 2012 -0800
files: compute_corrections.cxx constants.hxx main.cxx
description:
Compute dC_dv and test it
diff -r 5ad7ab8e65a2 -r 674bdfcc3cb8 compute_corrections.cxx
--- a/compute_corrections.cxx Wed Feb 22 13:30:45 2012 -0800
+++ b/compute_corrections.cxx Tue Feb 28 13:58:32 2012 -0800
@@ -68,7 +68,7 @@ void compute_corrections(const Model &mo
double x(i*h);
if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
{
- double Cx_new(0);
+ double Cx_new(0), dCx_dvx(0);
double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
@@ -89,10 +89,22 @@ void compute_corrections(const Model &mo
double dzx_xx_correction=
-(zx_jump + dx*dzx_x_jump
+ dx*dx*dzx_xx_jump/2)/(h*h);
- if(x>middle)
- dzx_xx_correction*=-1;
+ const int dzx_xx_sign(x>middle ? -1 : 1);
+ Cx_new=dzx_xx_sign*2*dzx_xx_correction;
- Cx_new=2*dzx_xx_correction;
+ double eta=exp(log_etax[i][j]);
+ double delta=(h-std::fabs(dx))/h;
+ double dvx_yy_dv=-delta/(h*h);
+ double dzy_yx_dv=-eta_jump*dvx_yy_dv;
+ double dvy_y_dv=-dzy_yx_dv*h/(min_eta+max_eta);
+ double dzx_x_dv=eta_jump*dvy_y_dv;
+ double dvx_dv=(4*eta*delta + 2*h*dzx_x_dv)/(3*(min_eta+max_eta));
+ double dp_x_dv=2*eta_jump*dvx_yy_dv;
+ double dzx_xx_dv=-eta_jump*dvx_yy_dv + dp_x_dv;
+ double dp_dv=-2*eta_jump*dvy_y_dv;
+
+ dCx_dvx=-dzx_xx_sign*2*(eta_jump*dvx_dv + dx*dzx_x_dv
+ + dx*dx*dzx_xx_dv/2)/(h*h);
if(x+h/2>middle && x-h/2<middle)
{
@@ -101,7 +113,33 @@ void compute_corrections(const Model &mo
Cx_new+=(p_jump + dx*dp_x_jump)/h;
Cx_new-=eta_jump*(dvy_y - dx*dvx_yy)/h;
+
+
+ dCx_dvx+= -(eta_jump*dvy_y_dv + dx*dzy_yx_dv)/h
+ + (dp_dv + dx*dp_x_dv)/h;
+
+ // if(j==Ny/2)
+ // std::cout << "Cx "
+ // << i << " "
+ // << j << " "
+ // << i*h << " "
+ // << eta << " "
+ // // << dzx_xx_sign << " "
+ // << zx[i][j] << " "
+ // << zy[i][j] << " "
+ // << Cx_new << " "
+ // << dCx_dvx/eta << " "
+ // // << dzx_xx_jump << " "
+ // // << dzx_xx_dv/eta << " "
+ // // << dzx_x_jump << " "
+ // // << dzx_x_dv/eta << " "
+ // // << zx_jump << " "
+ // // << eta_jump*dvx_dv/eta << " "
+ // << "\n";
+
+
}
+
update_correction(Cx_new,theta_correction,Cx[i][j]);
// Cx[i][j]+=theta_correction*(Cx_new-Cx[i][j]);
@@ -111,7 +149,7 @@ void compute_corrections(const Model &mo
x=(i+0.5)*h;
if(i<Nx && j!=0 && j!=Ny && x-h<middle && x+h>middle)
{
- double Cy_new(0);
+ double Cy_new(0), dCy_dvy(0);
double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
compute_v_on_interface(zx,zy,log_etax,
@@ -127,9 +165,23 @@ void compute_corrections(const Model &mo
double dzy_xx_jump=-3*eta_jump*dvy_yy;
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;
- Cy_new=dzy_xx_correction;
+ const int dzy_xx_sign(x>middle ? -1 : 1);
+ Cy_new=dzy_xx_sign*dzy_xx_correction;
+
+
+ double eta=exp(log_etay[i][j]);
+ double delta=(h-std::fabs(dx))/h;
+ double dvy_yy_dv=-delta/(h*h);
+ double dzx_yx_dv=-eta_jump*dvy_yy_dv;
+ double dvx_y_dv=(-h/(min_eta+max_eta))*dzx_yx_dv;
+ double dzy_x_dv=-eta_jump*dvx_y_dv;
+ double dvy_dv=
+ (4*eta*delta - 2*h*dzy_x_dv)/(3*(min_eta+max_eta));
+ double dp_y_dv=-2*eta_jump*dvy_yy_dv;
+ double dzy_xx_dv=-eta_jump*dvy_yy_dv + dp_y_dv;
+
+ dCy_dvy=-dzy_xx_sign*(eta_jump*dvy_dv - dx*dzy_x_dv
+ + dx*dx*dzy_xx_dv/2)/(h*h);
if(x+h/2>middle && x-h/2<middle)
{
@@ -138,6 +190,30 @@ void compute_corrections(const Model &mo
double dzx_yx_correction=
eta_jump*(dvx_y - dx*dvy_yy)/h;
Cy_new-=dzx_yx_correction;
+
+ dCy_dvy-=(eta_jump*dvx_y_dv + dx*dzx_yx_dv)/h;
+
+ if(j==Ny/2)
+ std::cout << "Cy "
+ << i << " "
+ << j << " "
+ << (i+0.5)*h << " "
+ << eta << " "
+ << dzy_xx_sign << " "
+ << zx[i][j] << " "
+ << zy[i][j] << " "
+ << Cy_new << " "
+ << dCy_dvy/eta << " "
+ // << dzy_xx_jump << " "
+ // << dzy_xx_dv/eta << " "
+ // << dzy_x_jump << " "
+ // << dzy_x_dv/eta << " "
+ // << zy_jump << " "
+ // << eta_jump*dvy_dv/eta << " "
+ << "\n";
+
+
+
}
update_correction(Cy_new,theta_correction,Cy[i][j]);
// Cy[i][j]+=theta_correction*(Cy_new-Cy[i][j]);
diff -r 5ad7ab8e65a2 -r 674bdfcc3cb8 constants.hxx
--- a/constants.hxx Wed Feb 22 13:30:45 2012 -0800
+++ b/constants.hxx Tue Feb 28 13:58:32 2012 -0800
@@ -1,14 +1,14 @@ const int N(128);
-const int N(128);
+const int N(64);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
-const double max_eta=1e1;
+const double max_eta=1e2;
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 + 0.00000001)/64 + 1/(2*Nx);
+const double middle=(25 + 0.00000001)/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);
diff -r 5ad7ab8e65a2 -r 674bdfcc3cb8 main.cxx
--- a/main.cxx Wed Feb 22 13:30:45 2012 -0800
+++ b/main.cxx Tue Feb 28 13:58:32 2012 -0800
@@ -34,8 +34,8 @@ int main()
/* Initial conditions */
- const int max_iterations=100000;
- int n_sweeps=10000000;
+ const int max_iterations=1;
+ int n_sweeps=1;
const double theta_mom=0.7;
const double theta_continuity=1.0;
const double theta_correction=0.01;
@@ -48,6 +48,16 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy,Cx,Cy,Cp);
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
+
+ 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;
+
+ zy[static_cast<int>(middle*Nx)][Ny/2]=1;
compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,1.0);
for(int iteration=0;iteration<max_iterations;++iteration)
@@ -492,6 +502,10 @@ int main()
if(sweep==0)
break;
}
+
+ write_vtk("zx",Nx+1,N,zx);
+ write_vtk("zy",Nx,N,zy);
+
for(int i=0;i<Nx+1;++i)
for(int j=0;j<Ny+1;++j)
{
More information about the CIG-COMMITS
mailing list