[cig-commits] commit: Compute all of dR_dz in a separate function.
Mercurial
hg at geodynamics.org
Thu Mar 29 18:35:08 PDT 2012
changeset: 135:40cd5b154dd2
user: Walter Landry <wlandry at caltech.edu>
date: Thu Mar 29 15:57:42 2012 -0700
files: compute_dC_dv.cxx compute_dR_dv.cxx main.cxx wscript
description:
Compute all of dR_dz in a separate function.
diff -r 546f3776347a -r 40cd5b154dd2 compute_dC_dv.cxx
--- a/compute_dC_dv.cxx Thu Mar 29 11:50:48 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-/* Computes the derivative of the correction w/respect to the
- augumented velocity (z). This is equal to the correction when
- z(i.j)=1 and everywhere z=0 and gravity=0, because everything is
- linear in z. */
-
-#include "constants.hxx"
-#include "compute_Cxyz.hxx"
-
-void compute_dC_dv(const double log_etax[Nx+1][Ny],
- const double log_etay[Nx][Ny+1],
- const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1],
- double dCx[Nx+1][Ny],
- double dCy[Nx][Ny+1])
-{
- double zx[Nx+1][Ny], zy[Nx][Ny+1];
- double fx[Nx+1][Ny], fy[Nx][Ny+1];
-
- for(int i=0;i<Nx+1;++i)
- for(int j=0;j<Ny;++j)
- zx[i][j]=fx[i][j]=0;
-
- for(int i=0;i<Nx;++i)
- for(int j=0;j<Ny+1;++j)
- zy[i][j]=fy[i][j]=0;
-
- for(int j=0;j<Ny;++j)
- for(int i=1;i<Nx;++i)
- {
- zx[i][j]=1;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 0,i,j,dCx[i][j]);
- zx[i][j]=0;
- }
-
- for(int j=1;j<Ny;++j)
- for(int i=0;i<Nx;++i)
- {
- zy[i][j]=1;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 1,i,j,dCy[i][j]);
- zy[i][j]=0;
- }
-
-}
-
-
diff -r 546f3776347a -r 40cd5b154dd2 compute_dR_dv.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_dR_dv.cxx Thu Mar 29 15:57:42 2012 -0700
@@ -0,0 +1,59 @@
+/* Computes the derivative of the residual and correction w/respect to the
+ augumented velocity (z). This is equal to the correction when
+ z(i.j)=1 and everywhere z=0 and gravity=0, because everything is
+ linear in z. */
+
+#include "constants.hxx"
+#include "compute_Cxyz.hxx"
+
+void compute_dR_dv(const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ double dRx[Nx+1][Ny],
+ double dRy[Nx][Ny+1])
+{
+ double zx[Nx+1][Ny], zy[Nx][Ny+1];
+ double fx[Nx+1][Ny], fy[Nx][Ny+1];
+
+ for(int i=0;i<Nx+1;++i)
+ for(int j=0;j<Ny;++j)
+ zx[i][j]=fx[i][j]=0;
+
+ for(int i=0;i<Nx;++i)
+ for(int j=0;j<Ny+1;++j)
+ zy[i][j]=fy[i][j]=0;
+
+ /* TODO: add in the terms due to the derivative of the viscosity */
+
+ for(int j=0;j<Ny;++j)
+ for(int i=1;i<Nx;++i)
+ {
+ zx[i][j]=1;
+ compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
+ 0,i,j,dRx[i][j]);
+ zx[i][j]=0;
+
+ if(j==0 || j==Ny-1)
+ dRx[i][j]-=5/(h*h);
+ else
+ dRx[i][j]-=6/(h*h);
+ }
+
+ for(int j=1;j<Ny;++j)
+ for(int i=0;i<Nx;++i)
+ {
+ zy[i][j]=1;
+ compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
+ 1,i,j,dRy[i][j]);
+ zy[i][j]=0;
+
+ if(i==0 || i==Nx-1)
+ dRy[i][j]-=5/(h*h);
+ else
+ dRy[i][j]-=6/(h*h);
+ }
+
+}
+
+
diff -r 546f3776347a -r 40cd5b154dd2 main.cxx
--- a/main.cxx Thu Mar 29 11:50:48 2012 -0700
+++ b/main.cxx Thu Mar 29 15:57:42 2012 -0700
@@ -21,12 +21,12 @@ extern void compute_Cp(const double zx[N
const int &i, const int &j,
double &Cp);
-extern void compute_dC_dv(const double log_etax[Nx+1][Ny],
+extern void compute_dR_dv(const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
- double dCx[Nx+1][Ny],
- double dCy[Nx][Ny+1]);
+ double dRx[Nx+1][Ny],
+ double dRy[Nx][Ny+1]);
int main()
@@ -39,25 +39,25 @@ int main()
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], fx[Nx+1][Ny], fy[Nx][Ny+1],
- distx[Nx+1][Ny], disty[Nx][Ny+1], dCx[Nx+1][Ny], dCy[Nx][Ny+1];
+ distx[Nx+1][Ny], disty[Nx][Ny+1], dRx[Nx+1][Ny], dRy[Nx][Ny+1];
double zx_new[Nx+1][Ny], zy_new[Nx][Ny+1], p_new[Nx][Ny];
- const bool check_initial(true);
+ const bool check_initial(false);
const bool jacobi(check_initial);
/* Initial conditions */
- const int n_sweeps(check_initial ? 1 : 1001);
+ const int n_sweeps(check_initial ? 1 : 10001);
const double theta_mom=0.7;
const double theta_continuity=1.0;
const double tolerance=1.0e-6;
- const int output_interval=10000000;
+ const int output_interval=1000;
Model model(Model::solCx);
initial(model,zx,zy,log_etax,log_etay,p,fx,fy,distx,disty);
- compute_dC_dv(log_etax,log_etay,distx,disty,dCx,dCy);
+ compute_dR_dv(log_etax,log_etay,distx,disty,dRx,dRy);
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
@@ -180,16 +180,11 @@ int main()
+ dlog_etaxy*zy_avg + dlog_etax*dzy_y
- dp_x - fx[i][j] + Cx;
- double dRx_dzx=-6/(h*h);
- if(j==0 || j==Ny-1)
- dRx_dzx=-5/(h*h);
- dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx[i][j];
-
Resid_x[i][j]=Rx;
if(!jacobi)
- zx[i][j]-=theta_mom*Rx/dRx_dzx;
+ zx[i][j]-=theta_mom*Rx/dRx[i][j];
else
- zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx_dzx;
+ zx_new[i][j]=zx[i][j]-theta_mom*Rx/dRx[i][j];
}
}
}
@@ -290,17 +285,12 @@ int main()
+ 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);
- if(i==0 || i==Nx-1)
- dRy_dzy=-5/(h*h);
- dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy[i][j];
-
+
Resid_y[i][j]=Ry;
if(!jacobi)
- zy[i][j]-=theta_mom*Ry/dRy_dzy;
+ zy[i][j]-=theta_mom*Ry/dRy[i][j];
else
- zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy_dzy;
+ zy_new[i][j]=zy[i][j]-theta_mom*Ry/dRy[i][j];
}
}
}
@@ -390,68 +380,28 @@ int main()
/* Fix vx */
if(i>0)
{
- double dlog_etaxx=
- ((log_etay[i][j+1] + log_etay[i][j])/2
- - 2*log_etax[i][j]
- + (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/(h*h/4);
-
- double dlog_etayy=
- ((log_etay[i][j+1] + log_etay[i-1][j+1])/2
- - 2*log_etax[i][j]
- + (log_etay[i][j] + log_etay[i-1][j])/2)/(h*h/4);
-
- if(model==Model::solCx)
- {
- dlog_etaxx=dlog_etayy=0;
- }
-
double Cx;
compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
0,i,j,Cx);
- double dRx_dzx=-6/(h*h);
- if(j==0 || j==Ny-1)
- dRx_dzx=-5/(h*h);
- dRx_dzx+=2*dlog_etaxx + dlog_etayy + dCx[i][j];
-
if(!jacobi)
- 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[i][j]);
else
- zx_new[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx_dzx);
+ zx_new[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx[i][j]);
max_x_resid=std::max(std::fabs(Resid_x[i][j]),max_x_resid);
}
/* Fix vy */
if(j>0)
{
- double dlog_etayy=
- ((log_etax[i+1][j] + log_etax[i][j])/2
- - 2*log_etay[i][j]
- + (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/(h*h/4);
-
- double dlog_etaxx=
- ((log_etax[i+1][j] + log_etax[i+1][j-1])/2
- - 2*log_etay[i][j]
- + (log_etax[i][j] + log_etax[i][j-1])/2)/(h*h/4);
-
- if(model==Model::solCx)
- {
- dlog_etaxx=dlog_etayy=0;
- }
-
double Cy;
compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
1,i,j,Cy);
- double dRy_dzy=-6/(h*h);
- if(i==0 || i==Nx-1)
- dRy_dzy=-5/(h*h);
- dRy_dzy+=2*dlog_etayy + dlog_etaxx + dCy[i][j];
-
if(!jacobi)
- 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[i][j]);
else
- zy_new[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy_dzy);
+ zy_new[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy[i][j]);
max_y_resid=std::max(fabs(Resid_y[i][j]),max_y_resid);
}
diff -r 546f3776347a -r 40cd5b154dd2 wscript
--- a/wscript Thu Mar 29 11:50:48 2012 -0700
+++ b/wscript Thu Mar 29 15:57:42 2012 -0700
@@ -12,7 +12,7 @@ def build(bld):
'compute_v_on_interface/compute_dv_dtt.cxx',
'compute_Cxyz.cxx',
'compute_Cp.cxx',
- 'compute_dC_dv.cxx',
+ 'compute_dR_dv.cxx',
'compute_jumps.cxx',
'Interface/Interface.cxx'],
More information about the CIG-COMMITS
mailing list