[cig-commits] commit: Compute the coefficients of the matrix beforehand. Speeds up by an
Mercurial
hg at geodynamics.org
Fri Mar 30 06:48:36 PDT 2012
changeset: 139:c63e54bcf2e9
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 30 06:48:19 2012 -0700
files: Coefficient.hxx compute_coefficients/compute_coefficients.cxx compute_coefficients/simplified_Rp.cxx compute_coefficients/simplified_Rx.cxx compute_coefficients/simplified_Ry.cxx compute_v_on_interface/compute_dv_dtt.cxx constants.hxx initial.cxx main.cxx wscript
description:
Compute the coefficients of the matrix beforehand. Speeds up by an
order of magnitude.
diff -r 4746e93fd167 -r c63e54bcf2e9 Coefficient.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Coefficient.hxx Fri Mar 30 06:48:19 2012 -0700
@@ -0,0 +1,13 @@
+#ifndef GAMR_COEFFICIENT_HXX
+#define GAMR_COEFFICIENT_HXX
+
+class Coefficient
+{
+public:
+ int i,j;
+ double val;
+ Coefficient(const int &I, const int &J, const double &V):
+ i(I), j(J), val(V) {}
+};
+
+#endif
diff -r 4746e93fd167 -r c63e54bcf2e9 compute_coefficients/compute_coefficients.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/compute_coefficients.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -0,0 +1,156 @@
+#include "../constants.hxx"
+#include "../Coefficient.hxx"
+#include <vector>
+#include <algorithm>
+#include "FTensor.hpp"
+
+template< typename T, class Allocator >
+void shrink_capacity(std::vector<T,Allocator>& v)
+{
+ std::vector<T,Allocator>(v.begin(),v.end()).swap(v);
+}
+
+double simplified_Rx(const int &i, const int &j,
+ 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 fx[Nx+1][Ny],
+ const double fy[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1]);
+
+double simplified_Ry(const int &i, const int &j,
+ 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 fx[Nx+1][Ny],
+ const double fy[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1]);
+
+double simplified_Rp(const int &i, const int &j,
+ 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 fx[Nx+1][Ny],
+ const double fy[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1]);
+
+void compute_coefficients(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],
+ std::vector<Coefficient> Cfx[2][Nx+1][Ny],
+ std::vector<Coefficient> Cfy[2][Nx][Ny+1],
+ std::vector<Coefficient> Cfp[2][Nx][Ny])
+{
+ double zx[Nx+1][Ny], zy[Nx][Ny+1], p[Nx][Ny];
+ double fx[Nx+1][Ny], fy[Nx][Ny+1];
+
+ for(auto &I: zx)
+ for(auto &J: I)
+ J=0;
+
+ for(auto &I: zy)
+ for(auto &J: I)
+ J=0;
+
+ for(auto &I: p)
+ for(auto &J: I)
+ J=0;
+
+ for(auto &I: fx)
+ for(auto &J: I)
+ J=0;
+
+ for(auto &I: fy)
+ for(auto &J: I)
+ J=0;
+
+ FTensor::Tensor2<int,2,2> max_ij(Nx,Ny-1,Nx-1,Ny);
+ for(int i=1;i<Nx;++i)
+ for(int j=0;j<Ny;++j)
+ {
+ for(int di=std::max(i-max_r,0);di<=std::min(i+max_r,max_ij(0,0));++di)
+ for(int dj=std::max(j-max_r,0);dj<=std::min(j+max_r,max_ij(0,1));++dj)
+ {
+ zx[di][dj]=1;
+ double R=simplified_Rx(i,j,zx,zy,log_etax,log_etay,fx,fy,
+ distx,disty);
+ if(R!=0)
+ Cfx[0][i][j].push_back(Coefficient(di,dj,R));
+ zx[di][dj]=0;
+ }
+ for(int di=std::max(i-max_r,0);di<=std::min(i+max_r,max_ij(1,0));++di)
+ for(int dj=std::max(j-max_r,0);dj<=std::min(j+max_r,max_ij(1,1));++dj)
+ {
+ zy[di][dj]=1;
+ double R=simplified_Rx(i,j,zx,zy,log_etax,log_etay,fx,fy,
+ distx,disty);
+ if(R!=0)
+ Cfx[1][i][j].push_back(Coefficient(di,dj,R));
+ zy[di][dj]=0;
+ }
+ shrink_capacity(Cfx[0][i][j]);
+ shrink_capacity(Cfx[1][i][j]);
+ }
+
+ for(int i=0;i<Nx;++i)
+ for(int j=1;j<Ny;++j)
+ {
+ for(int di=std::max(i-max_r,0);di<=std::min(i+max_r,max_ij(0,0));++di)
+ for(int dj=std::max(j-max_r,0);dj<=std::min(j+max_r,max_ij(0,1));++dj)
+ {
+ zx[di][dj]=1;
+ double R=simplified_Ry(i,j,zx,zy,log_etax,log_etay,fx,fy,
+ distx,disty);
+ if(R!=0)
+ Cfy[0][i][j].push_back(Coefficient(di,dj,R));
+ zx[di][dj]=0;
+ }
+ for(int di=std::max(i-max_r,0);di<=std::min(i+max_r,max_ij(1,0));++di)
+ for(int dj=std::max(j-max_r,0);dj<=std::min(j+max_r,max_ij(1,1));++dj)
+ {
+ zy[di][dj]=1;
+ double R=simplified_Ry(i,j,zx,zy,log_etax,log_etay,fx,fy,
+ distx,disty);
+ if(R!=0)
+ Cfy[1][i][j].push_back(Coefficient(di,dj,R));
+ zy[di][dj]=0;
+ }
+ shrink_capacity(Cfy[0][i][j]);
+ shrink_capacity(Cfy[1][i][j]);
+ }
+
+ for(int i=0;i<Nx;++i)
+ for(int j=0;j<Ny;++j)
+ {
+ for(int di=std::max(i-max_r,0);di<=std::min(i+max_r,max_ij(0,0));++di)
+ for(int dj=std::max(j-max_r,0);dj<=std::min(j+max_r,max_ij(0,1));++dj)
+ {
+ zx[di][dj]=1;
+ double R=simplified_Rp(i,j,zx,zy,log_etax,log_etay,fx,fy,
+ distx,disty);
+ if(R!=0)
+ Cfp[0][i][j].push_back(Coefficient(di,dj,R));
+ zx[di][dj]=0;
+ }
+ for(int di=std::max(i-max_r,0);di<=std::min(i+max_r,max_ij(1,0));++di)
+ for(int dj=std::max(j-max_r,0);dj<=std::min(j+max_r,max_ij(1,1));++dj)
+ {
+ zy[di][dj]=1;
+ double R=simplified_Rp(i,j,zx,zy,log_etax,log_etay,fx,fy,
+ distx,disty);
+ if(R!=0)
+ Cfp[1][i][j].push_back(Coefficient(di,dj,R));
+ zy[di][dj]=0;
+ }
+ shrink_capacity(Cfp[0][i][j]);
+ shrink_capacity(Cfp[1][i][j]);
+ }
+}
+
diff -r 4746e93fd167 -r c63e54bcf2e9 compute_coefficients/simplified_Rp.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/simplified_Rp.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -0,0 +1,30 @@
+/* Compute only the parts of Rp that depend on zx, zy */
+
+#include "../constants.hxx"
+
+extern void compute_Cp(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 distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const int &i, const int &j,
+ double &Cp);
+
+double simplified_Rp(const int &i, const int &j,
+ 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 fx[Nx+1][Ny],
+ const double fy[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1])
+{
+ double dzx_x=(zx[i+1][j] - zx[i][j])/h;
+ double dzy_y=(zy[i][j+1] - zy[i][j])/h;
+
+ double Cp;
+ compute_Cp(zx,zy,log_etax,log_etay,distx,disty,i,j,Cp);
+ return dzx_x + dzy_y + Cp;
+}
diff -r 4746e93fd167 -r c63e54bcf2e9 compute_coefficients/simplified_Rx.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/simplified_Rx.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -0,0 +1,40 @@
+/* Compute only the parts of Rx that depend on zx, zy */
+
+#include "../constants.hxx"
+#include "../compute_Cxyz.hxx"
+
+double simplified_Rx(const int &i, const int &j,
+ 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 fx[Nx+1][Ny],
+ const double fy[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1])
+{
+ double dzx_xx=
+ (zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
+
+ double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
+ - (zy[i][j] - zy[i-1][j]))/(h*h);
+
+ double dzx_yy;
+ if(j==0)
+ {
+ dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
+ }
+ else if(j==Ny-1)
+ {
+ dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
+ }
+ else
+ {
+ dzx_yy=(zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
+ }
+
+ double Cx;
+ compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,0,i,j,Cx);
+
+ return 2*dzx_xx + dzx_yy + dzy_xy + Cx;
+}
diff -r 4746e93fd167 -r c63e54bcf2e9 compute_coefficients/simplified_Ry.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_coefficients/simplified_Ry.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -0,0 +1,40 @@
+/* Compute only the parts of Ry that depend on zx, zy */
+
+#include "../constants.hxx"
+#include "../compute_Cxyz.hxx"
+
+double simplified_Ry(const int &i, const int &j,
+ 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 fx[Nx+1][Ny],
+ const double fy[Nx][Ny+1],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1])
+{
+ double dzy_yy=
+ (zy[i][j-1] - 2*zy[i][j] + zy[i][j+1])/(h*h);
+
+ double dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
+ - (zx[i][j] - zx[i][j-1]))/(h*h);
+
+ double dzy_xx;
+ if(i==0)
+ {
+ dzy_xx=(-zy[i][j] + zy[i+1][j])/(h*h);
+ }
+ else if(i==Nx-1)
+ {
+ dzy_xx=(zy[i-1][j] - zy[i][j])/(h*h);
+ }
+ else
+ {
+ dzy_xx=(zy[i-1][j] - 2*zy[i][j] + zy[i+1][j])/(h*h);
+ }
+
+ double Cy;
+ compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,1,i,j,Cy);
+
+ return 2*dzy_yy + dzy_xx + dzx_yx + Cy;
+}
diff -r 4746e93fd167 -r c63e54bcf2e9 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Thu Mar 29 20:11:29 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -25,8 +25,6 @@ void compute_dv_dtt(const double zx[Nx+1
FTensor::Tensor1<double,2> &dv,
FTensor::Tensor1<double,2> &v)
{
- const int max_r=4;
-
double length=h/std::max(std::abs(norm(0)),std::abs(norm(1)));
const FTensor::Index<'a',2> a;
diff -r 4746e93fd167 -r c63e54bcf2e9 constants.hxx
--- a/constants.hxx Thu Mar 29 20:11:29 2012 -0700
+++ b/constants.hxx Fri Mar 30 06:48:19 2012 -0700
@@ -15,6 +15,7 @@ const double middle=0.4;
const double middle=0.4;
const double h(1.0/N);
const double pi(atan(1.0)*4);
+const int max_r=4;
const double r2_inclusion=0.1;
const double epsilon=1;
diff -r 4746e93fd167 -r c63e54bcf2e9 initial.cxx
--- a/initial.cxx Thu Mar 29 20:11:29 2012 -0700
+++ b/initial.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -52,8 +52,9 @@ void initial(const Model &model, double
gsl_permutation_free(perm);
- const double theta(pi);
- // const double theta(0);
+ // const double theta(-0.5);
+ // const double theta(pi);
+ const double theta(0);
FTensor::Tensor2<double,2,2>
rot(std::cos(theta),std::sin(theta),-std::sin(theta),std::cos(theta));
FTensor::Tensor1<double,2> offset(0.5,0.5);
diff -r 4746e93fd167 -r c63e54bcf2e9 main.cxx
--- a/main.cxx Thu Mar 29 20:11:29 2012 -0700
+++ b/main.cxx Fri Mar 30 06:48:19 2012 -0700
@@ -6,6 +6,7 @@
#include "constants.hxx"
#include "write_vtk.hxx"
#include "compute_Cxyz.hxx"
+#include "Coefficient.hxx"
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],
@@ -28,6 +29,13 @@ extern void compute_dR_dv(const double l
double dRx[Nx+1][Ny],
double dRy[Nx][Ny+1]);
+void compute_coefficients(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],
+ std::vector<Coefficient> Cfx[2][Nx+1][Ny],
+ std::vector<Coefficient> Cfy[2][Nx][Ny+1],
+ std::vector<Coefficient> Cfp[2][Nx][Ny]);
int main()
{
@@ -41,6 +49,8 @@ int main()
p[Nx][Ny], dp[Nx][Ny], fx[Nx+1][Ny], fy[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];
+
+ std::vector<Coefficient> Cfx[2][Nx+1][Ny], Cfy[2][Nx][Ny+1], Cfp[2][Nx][Ny];
const bool check_initial(false);
const bool jacobi(check_initial);
@@ -58,6 +68,7 @@ int main()
initial(model,zx,zy,log_etax,log_etay,p,fx,fy,distx,disty);
compute_dR_dv(log_etax,log_etay,distx,disty,dRx,dRy);
+ compute_coefficients(log_etax,log_etay,distx,disty,Cfx,Cfy,Cfp);
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
@@ -97,88 +108,14 @@ int main()
{
if(i!=0 && i!=Nx)
{
- /* Do the finite difference parts */
-
- /* Derivatives of z */
- double dzx_xx=
- (zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
-
- double dzy_xy=((zy[i][j+1] - zy[i-1][j+1])
- - (zy[i][j] - zy[i-1][j]))/(h*h);
-
- double dzx_x=(zx[i+1][j] - zx[i-1][j])/(2*h);
- double dzy_y=(zy[i][j+1] - zy[i][j]
- + zy[i-1][j+1] - zy[i-1][j])/(2*h);
-
- double dzx_y, dzx_yy;
- if(j==0)
- {
- dzx_yy=(-zx[i][j] + zx[i][j+1])/(h*h);
- dzx_y=(zx[i][j+1] - zx[i][j])/(2*h);
- }
- else if(j==Ny-1)
- {
- dzx_yy=(zx[i][j-1] - zx[i][j])/(h*h);
- dzx_y=(zx[i][j] - zx[i][j-1])/(2*h);
- }
- else
- {
- dzx_yy=
- (zx[i][j-1] - 2*zx[i][j] + zx[i][j+1])/(h*h);
- dzx_y=(zx[i][j+1] - zx[i][j-1])/(2*h);
- }
-
- /* Derivatives of p and eta. */
-
double dp_x=(p[i][j]-p[i-1][j])/h;
- /* This plays a lot of games to reduce the
- size of the stencil. It may not be worth
- it, and it makes it more difficult to
- correct for the jumps. */
- 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_etax=
- ((log_etay[i][j+1] + log_etay[i][j])/2
- - (log_etay[i-1][j+1] + log_etay[i-1][j])/2)/h;
-
- 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);
-
- double dlog_etay=
- ((log_etay[i][j+1] + log_etay[i-1][j+1])
- - (log_etay[i][j] + log_etay[i-1][j]))/(2*h);
-
- double dlog_etaxy=
- ((log_etay[i][j+1] - log_etay[i-1][j+1])
- - (log_etay[i][j] - log_etay[i-1][j]))/(h*h);
-
- const double zy_avg=(zy[i][j] + zy[i-1][j] +
- zy[i][j+1] + zy[i-1][j+1])/4;
-
- /* Now do the jump conditions */
- if(model==Model::solCx)
- {
- dlog_etaxx=dlog_etax=dlog_etayy=dlog_etay=
- dlog_etaxy=0;
- }
-
- /* Compute the residual and update zx */
-
- double Cx;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 0,i,j,Cx);
-
- 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;
+ double Rx(0);
+ for(auto &c: Cfx[0][i][j])
+ Rx+=zx[c.i][c.j]*c.val;
+ for(auto &c: Cfx[1][i][j])
+ Rx+=zy[c.i][c.j]*c.val;
+ Rx-=dp_x + fx[i][j];
Resid_x[i][j]=Rx;
if(!jacobi)
@@ -207,85 +144,15 @@ int main()
{
if(j!=0 && j!=Ny)
{
- /* Do the finite difference parts */
-
- /* Derivatives of z */
- double dzy_yy=
- (zy[i][j-1] - 2*zy[i][j] + zy[i][j+1])/(h*h);
-
- double dzx_yx=((zx[i+1][j] - zx[i+1][j-1])
- - (zx[i][j] - zx[i][j-1]))/(h*h);
-
- double dzy_y=(zy[i][j+1] - zy[i][j-1])/(2*h);
- double dzx_x=(zx[i+1][j] - zx[i][j]
- + zx[i+1][j-1] - zx[i][j-1])/(2*h);
-
- double dzy_x, dzy_xx;
- if(i==0)
- {
- dzy_xx=(-zy[i][j] + zy[i+1][j])/(h*h);
- dzy_x=(zy[i+1][j] - zy[i][j])/(2*h);
- }
- else if(i==Nx-1)
- {
- dzy_xx=(zy[i-1][j] - zy[i][j])/(h*h);
- dzy_x=(zy[i][j] - zy[i-1][j])/(2*h);
- }
- else
- {
- dzy_xx=
- (zy[i-1][j] - 2*zy[i][j] + zy[i+1][j])/(h*h);
- dzy_x=(zy[i+1][j] - zy[i-1][j])/(2*h);
- }
-
- /* Derivatives of p and eta. */
-
double dp_y=(p[i][j]-p[i][j-1])/h;
- 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_etay=
- ((log_etax[i+1][j] + log_etax[i][j])/2
- - (log_etax[i+1][j-1] + log_etax[i][j-1])/2)/h;
+ double Ry(0);
+ for(auto &c: Cfy[0][i][j])
+ Ry+=zx[c.i][c.j]*c.val;
+ for(auto &c: Cfy[1][i][j])
+ Ry+=zy[c.i][c.j]*c.val;
+ Ry-=dp_y + fy[i][j];
- 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);
-
- double dlog_etax=
- ((log_etax[i+1][j] + log_etax[i+1][j-1])
- - (log_etax[i][j] + log_etax[i][j-1]))/(2*h);
-
- double dlog_etayx=
- ((log_etax[i+1][j] - log_etax[i+1][j-1])
- - (log_etax[i][j] - log_etax[i][j-1]))/(h*h);
-
- const double zx_avg=(zx[i][j] + zx[i][j-1] +
- zx[i+1][j] + zx[i+1][j-1])/4;
-
- /* Now do the jump conditions */
- if(model==Model::solCx)
- {
- dlog_etayy=dlog_etay=dlog_etaxx=dlog_etax=
- dlog_etayx=0;
- }
-
- /* Compute the residual and update zy */
-
- double Cy;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 1,i,j,Cy);
-
- 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;
-
Resid_y[i][j]=Ry;
if(!jacobi)
zy[i][j]-=theta_mom*Ry/dRy[i][j];
@@ -308,50 +175,40 @@ int main()
for(int j=0;j<Ny;++j)
for(int i=0;i<Nx;++i)
{
- double dzx_x=(zx[i+1][j] - zx[i][j])/h;
- double dzy_y=(zy[i][j+1] - zy[i][j])/h;
+ double Rc(0);
+ for(auto &c: Cfp[0][i][j])
+ Rc+=zx[c.i][c.j]*c.val;
+ for(auto &c: Cfp[1][i][j])
+ Rc+=zy[c.i][c.j]*c.val;
- double dlog_etax=(log_etax[i+1][j] - log_etax[i][j])/h;
- double dlog_etay=(log_etay[i][j+1] - log_etay[i][j])/h;
-
- double dlog_etaxx=
- (log_etax[i+1][j] - (log_etay[i][j+1] + log_etay[i][j])
- + log_etax[i][j])/(h*h/4);
- double dlog_etayy=
- (log_etay[i][j+1] - (log_etax[i+1][j] + log_etax[i][j])
- + log_etay[i][j])/(h*h/4);
-
- double zx_avg=(zx[i+1][j] + zx[i][j])/2;
- double zy_avg=(zy[i][j+1] + zy[i][j])/2;
-
- /* Apply the jump condition */
- if(model==Model::solCx)
- {
- dlog_etax=dlog_etay=dlog_etaxx=dlog_etayy=0;
- }
-
- double Cp;
- compute_Cp(zx,zy,log_etax,log_etay,distx,disty,i,j,Cp);
- double Rc=dzx_x + dzy_y + dlog_etax*zx_avg + dlog_etay*zy_avg
- + Cp;
-
double dRmp_dp=-1/h;
double dRmm_dp=1/h;
double dRc_dp(0);
if(i!=Nx-1)
- dRc_dp+=(1/h - dlog_etax/2)*dRmp_dp/dRx[i+1][j];
+ dRc_dp+=(1/h)*dRmp_dp/dRx[i+1][j];
if(i!=0)
- dRc_dp+=(-1/h - dlog_etax/2)*dRmm_dp/dRx[i][j];
+ dRc_dp+=(-1/h)*dRmm_dp/dRx[i][j];
if(j!=Ny-1)
- dRc_dp+=(1/h - dlog_etay/2)*dRmp_dp/dRy[i][j+1];
+ dRc_dp+=(1/h)*dRmp_dp/dRy[i][j+1];
if(j!=0)
- dRc_dp+=(-1/h - dlog_etay/2)*dRmm_dp/dRy[i][j];
+ dRc_dp+=(-1/h)*dRmm_dp/dRy[i][j];
- // double dRc_dp=4.0/6;
+ // if(i!=Nx-1)
+ // dRc_dp+=(1/h - dlog_etax/2)*dRmp_dp/dRx[i+1][j];
+
+ // if(i!=0)
+ // dRc_dp+=(-1/h - dlog_etax/2)*dRmm_dp/dRx[i][j];
+
+ // if(j!=Ny-1)
+ // dRc_dp+=(1/h - dlog_etay/2)*dRmp_dp/dRy[i][j+1];
+
+ // if(j!=0)
+ // dRc_dp+=(-1/h - dlog_etay/2)*dRmm_dp/dRy[i][j];
+
Resid_p[i][j]=Rc;
@@ -380,10 +237,6 @@ int main()
/* Fix vx */
if(i>0)
{
- double Cx;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 0,i,j,Cx);
-
if(!jacobi)
zx[i][j]+=(dp[i][j]-dp[i-1][j])/(h*dRx[i][j]);
else
@@ -394,10 +247,6 @@ int main()
/* Fix vy */
if(j>0)
{
- double Cy;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,
- 1,i,j,Cy);
-
if(!jacobi)
zy[i][j]+=(dp[i][j]-dp[i][j-1])/(h*dRy[i][j]);
else
diff -r 4746e93fd167 -r c63e54bcf2e9 wscript
--- a/wscript Thu Mar 29 20:11:29 2012 -0700
+++ b/wscript Fri Mar 30 06:48:19 2012 -0700
@@ -8,6 +8,10 @@ def build(bld):
features = ['cxx','cprogram'],
source = ['main.cxx',
'initial.cxx',
+ 'compute_coefficients/compute_coefficients.cxx',
+ 'compute_coefficients/simplified_Rx.cxx',
+ 'compute_coefficients/simplified_Ry.cxx',
+ 'compute_coefficients/simplified_Rp.cxx',
'compute_v_on_interface/compute_v_on_interface.cxx',
'compute_v_on_interface/compute_dv_dtt.cxx',
'compute_Cxyz.cxx',
More information about the CIG-COMMITS
mailing list