[cig-commits] commit: Rearrange so that Cp, Cx, and Cy are computed at each place during the iteration, rather than beforehand. Includes updated dC for higher order interpolation of v and testing stuff.
Mercurial
hg at geodynamics.org
Wed Feb 29 19:20:12 PST 2012
changeset: 63:bd123bf91338
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 29 15:30:17 2012 -0800
files: compute_Cp.cxx compute_Cx.cxx compute_Cy.cxx compute_corrections.cxx compute_v_on_interface.hxx constants.hxx initial.cxx main.cxx wscript
description:
Rearrange so that Cp, Cx, and Cy are computed at each place during the iteration, rather than beforehand. Includes updated dC for higher order interpolation of v and testing stuff.
diff -r c996563dfb3d -r bd123bf91338 compute_Cp.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_Cp.cxx Wed Feb 29 15:30:17 2012 -0800
@@ -0,0 +1,34 @@
+#include "constants.hxx"
+#include "Model.hxx"
+#include <iostream>
+#include "compute_v_on_interface.hxx"
+
+void compute_Cp(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cp)
+{
+ Cp=0;
+ if(model==Model::solCx)
+ {
+ const double x=(i+0.5)*h;
+ if(i<Nx && j<Ny && x+h/2>middle && x-h/2<middle)
+ {
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ compute_v_on_interface(zx,zy,log_etax,log_etay,
+ middle,j,0,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
+ double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
+ double zx_jump=eta_jump*vx;
+ double dzx_x_jump=eta_jump*dvy_y;
+
+ double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
+ Cp=-dzx_x_correction;
+ }
+ }
+}
diff -r c996563dfb3d -r bd123bf91338 compute_Cx.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_Cx.cxx Wed Feb 29 15:30:17 2012 -0800
@@ -0,0 +1,97 @@
+#include "constants.hxx"
+#include "Model.hxx"
+#include <iostream>
+#include "compute_v_on_interface.hxx"
+
+void compute_Cx(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cx, double &dCx_dzx)
+{
+ Cx=0;
+ dCx_dzx=0;
+ if(model==Model::solCx)
+ {
+ const double x(i*h);
+ if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
+ {
+ double dCx_dvx(0);
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ compute_v_on_interface(zx,zy,log_etax,
+ log_etay,
+ middle,j,0,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
+ double dx=(x>middle) ? (middle-(x-h))
+ : (middle-(x+h));
+ double zx_jump=eta_jump*vx;
+ double dzx_x_jump=eta_jump*dvy_y;
+ double p_jump=-2*dvy_y*eta_jump;
+ double dp_x_jump=2*dvx_yy*eta_jump;
+ double dzx_xx_jump=
+ -dvx_yy*eta_jump + dp_x_jump;
+ double dzx_xx_correction=
+ -(zx_jump + dx*dzx_x_jump
+ + dx*dx*dzx_xx_jump/2)/(h*h);
+ const int dzx_xx_sign(x>middle ? -1 : 1);
+ Cx=dzx_xx_sign*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 dvp_dv=delta/2 + delta*delta/2;
+ double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
+ double dvx_dv=(eta*(4*dvp_dv - dvpp_dv) + 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)
+ {
+ dx=(x>middle) ? ((x-h/2)-middle)
+ : ((x+h/2)-middle);
+
+ Cx+=(p_jump + dx*dp_x_jump)/h;
+ Cx-=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/8)
+ std::cout << "Cx "
+ << i << " "
+ << j << " "
+ << x << " "
+ << middle << " "
+ << zx[i][j] << " "
+ << Cx << " "
+ << dCx_dvx/eta << " "
+ << "\n";
+
+ }
+ if(j==Ny/8)
+ std::cout << "Cx out "
+ << i << " "
+ << j << " "
+ << x << " "
+ << middle << " "
+ << zx[i][j] << " "
+ << Cx << " "
+ << dCx_dvx/eta << " "
+ << "\n";
+
+ dCx_dzx=dCx_dvx/eta;
+ }
+ }
+}
diff -r c996563dfb3d -r bd123bf91338 compute_Cy.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_Cy.cxx Wed Feb 29 15:30:17 2012 -0800
@@ -0,0 +1,91 @@
+#include "constants.hxx"
+#include "Model.hxx"
+#include <iostream>
+#include "compute_v_on_interface.hxx"
+
+void compute_Cy(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cy, double &dCy_dzy)
+{
+ Cy=0;
+ dCy_dzy=0;
+ if(model==Model::solCx)
+ {
+ const double x=(i+0.5)*h;
+ if(i<Nx && j!=0 && j!=Ny && x-h<middle && x+h>middle)
+ {
+ double dCy_dvy(0);
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ compute_v_on_interface(zx,zy,log_etax,
+ log_etay,
+ middle,j,1,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
+ double dx=(x>middle) ? (middle-(x-h))
+ : (middle-(x+h));
+ double zy_jump=eta_jump*vy;
+ double dzy_x_jump=-eta_jump*dvx_y;
+ 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);
+ const int dzy_xx_sign(x>middle ? -1 : 1);
+ Cy=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 dvp_dv=delta/2 + delta*delta/2;
+ double dvpp_dv=-(1-delta)/2 + (1-delta)*(1-delta)/2;
+ double dvy_dv=
+ (eta*(4*dvp_dv - dvpp_dv) - 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)
+ {
+ dx=(x>middle) ? ((x-h/2)-middle)
+ : ((x+h/2)-middle);
+ double dzx_yx_correction=
+ eta_jump*(dvx_y - dx*dvy_yy)/h;
+ Cy-=dzx_yx_correction;
+
+ dCy_dvy-=(eta_jump*dvx_y_dv + dx*dzx_yx_dv)/h;
+
+ if(j==Ny/8)
+ std::cout << "Cy "
+ << i << " "
+ << j << " "
+ << x << " "
+ << middle << " "
+ << zy[i][j] << " "
+ << Cy << " "
+ << dCy_dvy/eta << " "
+ << "\n";
+ }
+
+ if(j==Ny/8)
+ std::cout << "Cy out "
+ << i << " "
+ << j << " "
+ << x << " "
+ << middle << " "
+ << zy[i][j] << " "
+ << Cy << " "
+ << dCy_dvy/eta << " "
+ << "\n";
+ dCy_dzy=dCy_dvy/eta;
+ }
+ }
+}
diff -r c996563dfb3d -r bd123bf91338 compute_corrections.cxx
--- a/compute_corrections.cxx Wed Feb 29 14:54:14 2012 -0800
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,245 +0,0 @@
-#include "constants.hxx"
-#include "Model.hxx"
-#include <iostream>
-
-extern void compute_v_on_interface(const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double &x, const int &j,
- const int &xyz,
- double &vx, double &vy,
- double &dvx_y, double &dvy_y,
- double &dvx_yy, double &dvy_yy);
-
-
-void update_correction(const double ¤t, const double &theta,
- double &old)
-{
- if(old==0)
- old=current;
- else
- {
- 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;
- }
-}
-
-void compute_corrections(const Model &model, const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- double Cx[Nx+1][Ny], double Cy[Nx][Ny+1],
- double Cp[Nx][Ny],
- const double &theta_correction)
-{
- if(model==Model::solCx)
- {
- for(int j=0;j<Ny+1;++j)
- for(int i=0;i<Nx+1;++i)
- {
- /* Cx */
- double x(i*h);
- if(j<Ny && i!=0 && i!=Nx && x-h<middle && x+h>middle)
- {
- double Cx_new(0), dCx_dvx(0);
-
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- compute_v_on_interface(zx,zy,log_etax,
- log_etay,
- middle,j,0,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- double dx=(x>middle) ? (middle-(x-h))
- : (middle-(x+h));
- double zx_jump=eta_jump*vx;
- double dzx_x_jump=eta_jump*dvy_y;
- double p_jump=-2*dvy_y*eta_jump;
- double dp_x_jump=2*dvx_yy*eta_jump;
- double dzx_xx_jump=
- -dvx_yy*eta_jump + dp_x_jump;
- double dzx_xx_correction=
- -(zx_jump + dx*dzx_x_jump
- + dx*dx*dzx_xx_jump/2)/(h*h);
- const int dzx_xx_sign(x>middle ? -1 : 1);
- Cx_new=dzx_xx_sign*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)
- {
- dx=(x>middle) ? ((x-h/2)-middle)
- : ((x+h/2)-middle);
-
- 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]);
- }
-
- /* Cy */
- x=(i+0.5)*h;
- if(i<Nx && j!=0 && j!=Ny && x-h<middle && x+h>middle)
- {
- 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,
- log_etay,
- middle,j,1,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- double dx=(x>middle) ? (middle-(x-h))
- : (middle-(x+h));
- double zy_jump=eta_jump*vy;
- double dzy_x_jump=-eta_jump*dvx_y;
- 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);
- 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)
- {
- dx=(x>middle) ? ((x-h/2)-middle)
- : ((x+h/2)-middle);
- 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]);
- }
-
- /* Cp */
- if(i<Nx && j<Ny && x+h/2>middle && x-h/2<middle)
- {
- double Cp_new(0);
- double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
-
- compute_v_on_interface(zx,zy,log_etax,log_etay,
- middle,j,0,
- vx,vy,dvx_y,dvy_y,
- dvx_yy,dvy_yy);
-
- double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
- double zx_jump=eta_jump*vx;
- double dzx_x_jump=eta_jump*dvy_y;
-
- double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
- Cp_new=-dzx_x_correction;
-
- update_correction(Cp_new,theta_correction,Cp[i][j]);
- // Cp[i][j]+=theta_correction*(Cp_new-Cp[i][j]);
- }
- }
- }
-}
diff -r c996563dfb3d -r bd123bf91338 compute_v_on_interface.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface.hxx Wed Feb 29 15:30:17 2012 -0800
@@ -0,0 +1,14 @@
+#ifndef GAMR_COMPUTE_V_ON_INTERFACE_HXX
+#define GAMR_COMPUTE_V_ON_INTERFACE_HXX
+
+extern void compute_v_on_interface(const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double &x, const int &j,
+ const int &xyz,
+ double &vx, double &vy,
+ double &dvx_y, double &dvy_y,
+ double &dvx_yy, double &dvy_yy);
+
+#endif
diff -r c996563dfb3d -r bd123bf91338 constants.hxx
--- a/constants.hxx Wed Feb 29 14:54:14 2012 -0800
+++ b/constants.hxx Wed Feb 29 15:30:17 2012 -0800
@@ -1,4 +1,4 @@ const int N(64);
-const int N(64);
+const int N(1024);
const int Nx(N);
const int Ny(N);
const double min_eta=1;
@@ -7,8 +7,8 @@ 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 c996563dfb3d -r bd123bf91338 initial.cxx
--- a/initial.cxx Wed Feb 29 14:54:14 2012 -0800
+++ b/initial.cxx Wed Feb 29 15:30:17 2012 -0800
@@ -8,8 +8,7 @@
void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
- double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
- double Cx[Nx+1][Ny], double Cy[Nx][Ny+1], double Cp[Nx][Ny])
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1])
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
@@ -55,7 +54,6 @@ void initial(const Model &model, double
{
if(i<Nx && j<Ny)
{
- Cp[i][j]=0;
double x((i+0.5)*h), y((j+0.5)*h), r2(x*x+y*y);
switch(model)
{
@@ -104,7 +102,6 @@ void initial(const Model &model, double
}
if(j<Ny)
{
- Cx[i][j]=0;
double x(i*h), y((j+0.5)*h), r2(x*x+y*y);
zx[i][j]=0;
fx[i][j]=0;
@@ -181,7 +178,6 @@ void initial(const Model &model, double
if(i<Nx)
{
- Cy[i][j]=0;
double x((i+0.5)*h), y(j*h), r2(x*x+y*y);
zy[i][j]=0;
diff -r c996563dfb3d -r bd123bf91338 main.cxx
--- a/main.cxx Wed Feb 29 14:54:14 2012 -0800
+++ b/main.cxx Wed Feb 29 15:30:17 2012 -0800
@@ -8,17 +8,29 @@
extern void initial(const Model &model, double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
- double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
- double Cx[Nx+1][Ny], double Cy[Nx][Ny+1],
- double Cp[Nx][Ny]);
+ double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1]);
-extern void compute_corrections(const Model &model, const double zx[Nx+1][Ny],
- const double zy[Nx][Ny+1],
- const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- double Cx[Nx+1][Ny], double Cy[Nx][Ny+1],
- double Cp[Nx][Ny],
- const double &theta_correction);
+extern void compute_Cx(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cx, double &dCx_dzx);
+
+extern void compute_Cy(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cy, double &dCy_dzy);
+
+extern void compute_Cp(const Model &model, const double zx[Nx+1][Ny],
+ const double zy[Nx][Ny+1],
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cp);
+
int main()
{
@@ -29,8 +41,7 @@ int main()
with z, it is easy to compute v=z/eta. */
double zx[Nx+1][Ny], zy[Nx][Ny+1], log_etax[Nx+1][Ny], log_etay[Nx][Ny+1],
- p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1],
- Cx[Nx+1][Ny], Cy[Nx][Ny+1], Cp[Nx][Ny];
+ p[Nx][Ny], dp[Nx][Ny], div[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1];
/* Initial conditions */
@@ -38,16 +49,16 @@ int main()
int n_sweeps=1;
const double theta_mom=0.7;
const double theta_continuity=1.0;
- const double theta_correction=0.01;
const double tolerance=1.0e-6;
const int output_interval=100000000;
Model model(Model::solCx);
- initial(model,zx,zy,log_etax,log_etay,p,fx,fy,Cx,Cy,Cp);
+ 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];
+
for(int i=0;i<Nx+1;++i)
for(int j=0;j<Ny;++j)
@@ -57,22 +68,18 @@ int main()
for(int j=0;j<Ny+1;++j)
zy[i][j]=0;
- zy[static_cast<int>(middle*Nx)][Ny/2]=1;
+ // zx[static_cast<int>(middle/h)][Ny/8]=1;
+ zy[static_cast<int>(middle/h)-1][Ny/8]=1;
- compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,1.0);
for(int iteration=0;iteration<max_iterations;++iteration)
{
- compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
/* Smoothing sweeps */
// if(iteration==max_iterations-1)
// n_sweeps=1;
- double max_x_resid_previous(0), max_y_resid_previous(0), max_p_resid_previous(0);
int sweep;
for(sweep=0;sweep<n_sweeps;++sweep)
{
- // compute_corrections(model,zx,zy,log_etax,log_etay,Cx,Cy,Cp,theta_correction);
-
/* zx */
for(int rb=0;rb<2;++rb)
{
@@ -156,37 +163,19 @@ int main()
/* Compute the residual and update zx */
-
- // if(j==Ny/8 && i*h>middle-3*h && i*h<middle+3*h)
- // std::cout << "Rx "
- // << i << " "
- // << j << " "
- // << i*h << " "
- // << middle << " "
- // << dzx_xx << " "
- // << Cx[i][j] << " "
- // << (2*dzx_xx + dzx_yy + dzy_xy
- // + 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
- // + dlog_etayy*zx[i][j] + dlog_etay*dzx_y
- // + dlog_etaxy*zy_avg + dlog_etax*dzy_y
- // - dp_x - fx[i][j]) << " "
- // << (2*dzx_xx + dzx_yy + dzy_xy
- // + 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
- // + dlog_etayy*zx[i][j] + dlog_etay*dzx_y
- // + dlog_etaxy*zy_avg + dlog_etax*dzy_y
- // - dp_x - fx[i][j] + Cx[i][j]) << " "
- // << "\n";
+ double Cx, dCx_dzx;
+ compute_Cx(model,zx,zy,log_etax,log_etay,i,j,Cx,dCx_dzx);
double Rx=2*dzx_xx + dzx_yy + dzy_xy
+ 2*(dlog_etaxx*zx[i][j] + dlog_etax*dzx_x)
+ dlog_etayy*zx[i][j] + dlog_etay*dzx_y
+ dlog_etaxy*zy_avg + dlog_etax*dzy_y
- - dp_x - fx[i][j] + Cx[i][j];
+ - dp_x - fx[i][j] + Cx;
double dRx_dzx=-6/(h*h) + 2*dlog_etaxx + dlog_etayy;
Resid_x[i][j]=Rx;
- zx[i][j]-=theta_mom*Rx/dRx_dzx;
+ // zx[i][j]-=theta_mom*Rx/dRx_dzx;
}
}
}
@@ -197,12 +186,12 @@ int main()
std::stringstream ss;
ss << "zx_resid" << sweep;
// 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);
- }
+ // if(sweep==0)
+ // {
+ // ss.str("");
+ // ss << "Cx" << iteration;
+ // write_vtk(ss.str(),Nx+1,N,Cx);
+ // }
}
/* zy */
@@ -282,29 +271,40 @@ int main()
dlog_etayx=0;
}
- // if(j==Ny/8 && i*h>middle-3*h && i*h<middle+3*h)
+ /* Compute the residual and update zy */
+
+ double Cy, dCy_dzy;
+ compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
+ dCy_dzy);
+ double Ry=2*dzy_yy + dzy_xx + dzx_yx
+ + 2*(dlog_etayy*zy[i][j] + dlog_etay*dzy_y)
+ + dlog_etaxx*zy[i][j] + dlog_etax*dzy_x
+ + dlog_etayx*zx_avg + dlog_etay*dzx_x
+ - dp_y - fy[i][j] + Cy;
+
+ double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
+
+ // if(j==Ny/8 && i*h>middle-10*h && i*h<middle+10*h)
// std::cout << "Ry "
// << i << " "
// << j << " "
// << i*h+h/2 << " "
// << middle << " "
- // << dzy_xx << " "
- // << Cy[i][j] << " "
- // << dzy_xx+Cy[i][j] << " "
+
+ // // << 2*dzy_yy << " "
+ // // << dzy_xx << " "
+ // // << dzx_yx << " "
+ // // << Ry << " "
+ // // << Cy << " "
+
+ // << dzy_xx + Cy << " "
+
// << "\n";
- /* Compute the residual and update zy */
- double Ry=2*dzy_yy + dzy_xx + dzx_yx
- + 2*(dlog_etayy*zy[i][j] + dlog_etay*dzy_y)
- + dlog_etaxx*zy[i][j] + dlog_etax*dzy_x
- + dlog_etayx*zx_avg + dlog_etay*dzx_x
- - dp_y - fy[i][j] + Cy[i][j];
-
- double dRy_dzy=-6/(h*h) + 2*dlog_etayy + dlog_etaxx;
Resid_y[i][j]=Ry;
- zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ // zy[i][j]-=theta_mom*Ry/dRy_dzy;
}
}
}
@@ -315,12 +315,12 @@ int main()
std::stringstream ss;
ss << "zy_resid" << sweep;
// write_vtk(ss.str(),Nx,N,Resid_y);
- if(sweep==0)
- {
- ss.str("");
- ss << "Cy" << iteration;
- write_vtk(ss.str(),Nx,N,Cy);
- }
+ // if(sweep==0)
+ // {
+ // ss.str("");
+ // ss << "Cy" << iteration;
+ // write_vtk(ss.str(),Nx,N,Cy);
+ // }
}
/* Pressure update */
@@ -350,7 +350,10 @@ int main()
dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
}
- double Rc=dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg + Cp[i][j];
+ double Cp;
+ compute_Cp(model,zx,zy,log_etax,log_etay,i,j,Cp);
+ double Rc=
+ dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg + Cp;
double dRc_dzxp=1/h - dlog_etax/2;
double dRc_dzxm=-1/h - dlog_etax/2;
@@ -370,7 +373,7 @@ int main()
Resid_p[i][j]=Rc;
dp[i][j]=-theta_continuity*Rc/dRc_dp;
- p[i][j]+=dp[i][j];
+ // p[i][j]+=dp[i][j];
}
if(sweep%output_interval==0)
@@ -378,12 +381,12 @@ int main()
std::stringstream ss;
ss << "p_resid" << sweep;
// write_vtk(ss.str(),Nx,N,Resid_p);
- if(sweep==0)
- {
- ss.str("");
- ss << "Cp" << iteration;
- write_vtk(ss.str(),Nx,N,Cp);
- }
+ // if(sweep==0)
+ // {
+ // ss.str("");
+ // ss << "Cp" << iteration;
+ // write_vtk(ss.str(),Nx,N,Cp);
+ // }
}
/* Velocity Fix */
@@ -413,7 +416,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);
}
@@ -437,7 +440,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);
}
@@ -446,7 +449,8 @@ int main()
if(sweep%1000==0)
{
- std::cout << iteration << " "
+ std::cout << "sweep "
+ << iteration << " "
<< sweep << " "
<< max_x_resid << " "
<< max_y_resid << " "
@@ -495,9 +499,6 @@ int main()
<< "\n";
break;
}
- max_x_resid_previous=max_x_resid;
- max_y_resid_previous=max_y_resid;
- max_p_resid_previous=max_p_resid;
}
if(sweep==0)
break;
diff -r c996563dfb3d -r bd123bf91338 wscript
--- a/wscript Wed Feb 29 14:54:14 2012 -0800
+++ b/wscript Wed Feb 29 15:30:17 2012 -0800
@@ -12,7 +12,9 @@ def build(bld):
source = ['main.cxx',
'initial.cxx',
'compute_v_on_interface.cxx',
- 'compute_corrections.cxx'],
+ 'compute_Cx.cxx',
+ 'compute_Cy.cxx',
+ 'compute_Cp.cxx'],
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-D_GLIBCXX_DEBUG','-Wall'],
More information about the CIG-COMMITS
mailing list