[cig-commits] commit: Compute Cy using compute_Cxyz
Mercurial
hg at geodynamics.org
Mon Mar 12 18:01:10 PDT 2012
changeset: 85:e4b405469c16
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Mon Mar 12 18:00:55 2012 -0700
files: compute_Cxyz.cxx compute_Cy.cxx compute_v_on_interface.cxx main.cxx wscript
description:
Compute Cy using compute_Cxyz
diff -r 14f98b7e3c8c -r e4b405469c16 compute_Cxyz.cxx
--- a/compute_Cxyz.cxx Mon Mar 12 14:32:21 2012 -0700
+++ b/compute_Cxyz.cxx Mon Mar 12 18:00:55 2012 -0700
@@ -38,7 +38,7 @@ void compute_Cxyz(const Model &model, co
FTensor::Tensor1<int,2> xyz(0,0), dir(0,0);
xyz(d)=1;
dir(dd)=1;
- compute_v_on_interface(zx,zy,log_etax,log_etay,interface.pos,i,j,0,
+ compute_v_on_interface(zx,zy,log_etax,log_etay,interface.pos,i,j,d,
v,dv,ddv,v_dv,dv_dv,ddv_dv);
FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
@@ -69,8 +69,10 @@ void compute_Cxyz(const Model &model, co
double dz_dd_correction=
-(z_jump - dx*dz_d_jump + dx*dx*ddz_dd_jump/2)/(h*h);
- const int dz_dd_factor((pos(dd)>interface.pos(dd) ? -1 : 1)*(d==dd ? 2 : 1));
- C=dz_dd_factor*dz_dd_correction;
+ const int dz_dd_factor((pos(dd)>interface.pos(dd) ? -1 : 1)
+ *(d==dd ? 2 : 1));
+
+ C+=dz_dd_factor*dz_dd_correction;
FTensor::Tensor1<double,2> dp_dv(2*eta_jump*ddv_dv(n),
-2*eta_jump*ddv_dv(t));
@@ -85,7 +87,7 @@ void compute_Cxyz(const Model &model, co
ddz_dv(a,n,n)=-ddz_dv(a,t,t) + dp_dv(a);
ddz_dv(a,n,t)=-eta_jump*ddv_dv(b)*ddel(a,b);
- dC_dv=-dz_dd_factor
+ dC_dv-=dz_dd_factor
*(eta_jump*v_dv(a)*xyz(a)
- dx*dz_dv(a,b)*xyz(a)*dir(b)
+ dx*dx*ddz_dv(a,b,c)*xyz(a)*dir(b)*dir(c)/2)/(h*h);
@@ -119,10 +121,10 @@ void compute_Cxyz(const Model &model, co
switch(d)
{
case 0:
- dC_dz=dC_dv*exp(-log_etax[i][j]);
+ dC_dz+=dC_dv*exp(-log_etax[i][j]);
break;
case 1:
- dC_dz=dC_dv*exp(-log_etay[i][j]);
+ dC_dz+=dC_dv*exp(-log_etay[i][j]);
break;
default:
abort();
diff -r 14f98b7e3c8c -r e4b405469c16 compute_Cy.cxx
--- a/compute_Cy.cxx Mon Mar 12 14:32:21 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-#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,i,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=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
- 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;
- }
- dCy_dzy=dCy_dvy/eta;
- }
- }
-}
diff -r 14f98b7e3c8c -r e4b405469c16 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Mon Mar 12 14:32:21 2012 -0700
+++ b/compute_v_on_interface.cxx Mon Mar 12 18:00:55 2012 -0700
@@ -284,8 +284,39 @@ void compute_v_on_interface(const double
dv(0)=compute_dv_y(zx,log_etax,dzx_yx_jump(ddv(1)),dx,ix,j-1);
v(1)=compute_v(zy,log_etay,dzy_x_jump(dv(0)),dy,iy,j);
- v(0)=std::numeric_limits<double>::infinity();
- dv(1)=std::numeric_limits<double>::infinity();
- ddv(0)=std::numeric_limits<double>::infinity();
+ v(0)=0;
+ dv(1)=0;
+ ddv(0)=0;
+
+ double eta=exp(log_etay[i][j]);
+ double x_i((i+0.5)*h);
+ double delta_x=(x_i>middle) ? (middle-(x_i-h))
+ : (middle-(x_i+h));
+ double delta=(h-std::fabs(delta_x))/h;
+ double dvy_yy_dv=-2*eta*eta/((h*h)*(min_eta*min_eta + max_eta*max_eta));
+ 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));
+
+ ddv_dv(0)=0;
+ ddv_dv(1)=dvy_yy_dv;
+ dv_dv(0)=dvx_y_dv;
+ dv_dv(1)=0;
+ v_dv(0)=0;
+ v_dv(1)=dvy_dv;
+
+
+ // if(j==1)
+ // std::cout << "v "
+ // << ddv_dv(1) << " "
+ // << dv_dv(0) << " "
+ // << v_dv(1) << " "
+ // << delta << " "
+ // << dy << " "
+ // << "\n";
}
}
diff -r 14f98b7e3c8c -r e4b405469c16 main.cxx
--- a/main.cxx Mon Mar 12 14:32:21 2012 -0700
+++ b/main.cxx Mon Mar 12 18:00:55 2012 -0700
@@ -18,13 +18,6 @@ extern void compute_Cxyz(const Model &mo
const int &i, const int &j,
double &C, double &dCx_dz);
-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],
@@ -48,7 +41,7 @@ int main()
const bool jacobi(false);
/* Initial conditions */
- int n_sweeps=100000000;
+ int n_sweeps=1000000;
const double theta_mom=0.7;
const double theta_continuity=1.0;
const double tolerance=1.0e-6;
@@ -272,8 +265,8 @@ int main()
/* Compute the residual and update zy */
double Cy, dCy_dzy;
- compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
- dCy_dzy);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,1,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
@@ -425,8 +418,8 @@ int main()
}
double Cy, dCy_dzy;
- compute_Cy(model,zx,zy,log_etax,log_etay,i,j,Cy,
- dCy_dzy);
+ compute_Cxyz(model,zx,zy,log_etax,log_etay,1,i,j,Cy,dCy_dzy);
+
double dRy_dzy=-6/(h*h);
if(i==0 || i==Nx-1)
dRy_dzy=-5/(h*h);
diff -r 14f98b7e3c8c -r e4b405469c16 wscript
--- a/wscript Mon Mar 12 14:32:21 2012 -0700
+++ b/wscript Mon Mar 12 18:00:55 2012 -0700
@@ -10,7 +10,6 @@ def build(bld):
'initial.cxx',
'compute_v_on_interface.cxx',
'compute_Cxyz.cxx',
- 'compute_Cy.cxx',
'compute_Cp.cxx'],
target = 'Gamr_disc',
More information about the CIG-COMMITS
mailing list