[cig-commits] commit: Compute interfaces in a more general way suitable for curved interfaces.
Mercurial
hg at geodynamics.org
Fri Apr 13 16:41:18 PDT 2012
changeset: 152:8945333f6754
user: Walter Landry <wlandry at caltech.edu>
date: Fri Apr 13 15:46:54 2012 -0700
files: compute_coefficients/Interface.hxx compute_coefficients/Interface/Interface.cxx compute_coefficients/compute_Cp.cxx compute_coefficients/compute_Cxyz.cxx compute_coefficients/compute_Cxyz.hxx compute_coefficients/compute_coefficients.cxx compute_coefficients/simplified_Rp.cxx compute_coefficients/simplified_Rx.cxx compute_coefficients/simplified_Ry.cxx initial.cxx main.cxx wscript
description:
Compute interfaces in a more general way suitable for curved interfaces.
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/Interface.hxx
--- a/compute_coefficients/Interface.hxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/Interface.hxx Fri Apr 13 15:46:54 2012 -0700
@@ -5,51 +5,240 @@
#include "sign.hxx"
#include "FTensor.hpp"
#include <iostream>
+#include <cassert>
class Interface
{
public:
- FTensor::Tensor1<double,2> grid_pos, dd_pos[2], corner_pos[2][2],
- anticorner_pos[2], dd_norm[2], corner_norm[2][2], anticorner_norm[2],
- dd_tangent[2], corner_tangent[2][2], anticorner_tangent[2];
- bool intersect_dd[2], intersect_sides[2], intersect_corner[2][2],
- intersect_anticorner[2][2];
- bool intersect_dp;
+ FTensor::Tensor1<double,2> grid_pos,
+ dd_pos[2][2], dd_norm[2][2], dd_tangent[2][2],
+ sides_pos[2][2], sides_norm[2][2], sides_tangent[2][2],
+ corner_pos[2][2], corner_norm[2][2], corner_tangent[2][2],
+ anticorner_pos[2], anticorner_norm[2], anticorner_tangent[2],
+ dp_pos[2], dp_norm[2], dp_tangent[2];
+ bool intersect_dd[2][2], intersect_sides[2][2], intersect_corner[2][2],
+ intersect_anticorner[2][2], intersect_dp[2];
+
+ double eta_jump_dd[2][2], eta_jump_sides[2][2], eta_jump_corner[2][2],
+ eta_jump_anticorner[2][2], eta_jump_dp[2];
Interface(const int &d, const int &i, const int &j,
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1]);
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
+ const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny], const double dist_edge[Nx+1][Ny+1]);
bool intersects() const
{
- return intersect_dd[0] || intersect_dd[1]
- || intersect_sides[0] || intersect_sides[1] || intersect_dp
+ return intersect_dd[0][0] || intersect_dd[0][1]
+ || intersect_dd[1][0] || intersect_dd[1][1]
+ || intersect_sides[0][0] || intersect_sides[0][1]
+ || intersect_sides[1][0] || intersect_sides[1][1]
|| intersect_corner[0][0] || intersect_corner[0][1]
- || intersect_corner[1][0] || intersect_corner[1][1];
+ || intersect_corner[1][0] || intersect_corner[1][1]
+ || intersect_anticorner[0][0] || intersect_anticorner[0][1]
+ || intersect_anticorner[1][0] || intersect_anticorner[1][1]
+ || intersect_dp[0] || intersect_dp[1];
}
-
- template <int ny> void compute_xy(const int &i, const int &j,
- const double ¢er_dist,
- const double dist[][ny])
+
+ static double intersection(const double &u_m, const double &u_0,
+ const double &u_p, const double &H)
{
+ /* This assumes that there is only one intersection between u_p and u_m */
+
+ double du=(u_p-u_m)/2;
+ double ddu=u_p-2*u_0+u_m;
+
+ if(std::fabs(ddu)>1e-8*(std::fabs(u_p)+std::fabs(u_0)+std::fabs(u_p)))
+ {
+ assert(du*du-2*u_0*ddu>=0);
+ double tmp=sqrt(du*du-2*u_0*ddu);
+
+ double p=(-du+tmp)/ddu;
+ double m=(-du-tmp)/ddu;
+
+ if(m>-1)
+ return m*H;
+ return p*H;
+ }
+ else
+ {
+ return -H*u_0/du;
+ }
+ }
+
+ template <int ny> void compute_dd(const int &d, const int &i, const int &j,
+ const double log_eta[][ny],
+ const double dist[][ny],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
+ {
+ const int nx=Nx+1-d;
+
+ /* TODO: fix if we have Neumann boundaries. Also, this may cut
+ off an interface at the boundary, because it assumes that if you
+ are at i==0, then the interface does not cut below the i=0
+ plane. */
+ intersect_dd[0][0]=(i!=0 && sign(dist[i][j])!=sign(dist[i-1][j]));
+ intersect_dd[0][1]=(i!=nx-1 && sign(dist[i][j])!=sign(dist[i+1][j]));
+
+ intersect_dd[1][0]=(j!=0 && sign(dist[i][j])!=sign(dist[i][j-1]));
+ intersect_dd[1][1]=(j!=ny-1 && sign(dist[i][j])!=sign(dist[i][j+1]));
+
+ /* TODO: Fix for variable viscosity */
+
+ if(intersect_dd[0][0])
+ eta_jump_dd[0][0]=exp(log_eta[i][j])-exp(log_eta[i-1][j]);
+ if(intersect_dd[0][1])
+ eta_jump_dd[0][1]=exp(log_eta[i][j])-exp(log_eta[i+1][j]);
+
+ if(intersect_dd[1][0])
+ eta_jump_dd[1][0]=exp(log_eta[i][j])-exp(log_eta[i][j-1]);
+ if(intersect_dd[1][1])
+ eta_jump_dd[1][1]=exp(log_eta[i][j])-exp(log_eta[i][j+1]);
+
+ for(int pm=0;pm<2;++pm)
+ {
+ for(int dd=0;dd<2;++dd)
+ {
+ if(intersect_dd[dd][pm])
+ {
+ dd_pos[dd][pm]((dd+1)%2)=grid_pos((dd+1)%2);
+ if(dd==d)
+ {
+ dd_pos[dd][pm](dd)=grid_pos(dd)
+ + (2*pm-1)*(intersection
+ (dist[i][j],
+ dist_cell[i+(pm-1)*(d==0)][j+(pm-1)*(d==1)],
+ dist[i+(2*pm-1)*(dd==0)][j+(2*pm-1)*(dd==1)],
+ h/2)
+ + h/2);
+ }
+ else
+ {
+ dd_pos[dd][pm](dd)=grid_pos(dd)
+ + (2*pm-1)*(intersection
+ (dist[i][j],
+ dist_edge[i+pm*(dd==0)][j+pm*(dd==1)],
+ dist[i+(2*pm-1)*(dd==0)][j+(2*pm-1)*(dd==1)],
+ h/2)
+ +h/2);
+ }
+ }
+ }
+ }
+ }
+
+ template <int ny> void compute_dp(const int &d, const int &i, const int &j,
+ const double log_eta[][ny],
+ const double eta_cell[Nx][Ny],
+ const double dist[][ny],
+ const double dist_cell[Nx][Ny])
+ {
+ for(int pm=0;pm<2;++pm)
+ {
+ int ip=i+(pm-1)*(d==0);
+ int jp=j+(pm-1)*(d==1);
+
+ intersect_dp[pm]=(sign(dist_cell[ip][jp])!=sign(dist[i][j]));
+
+ if(intersect_dp[pm])
+ {
+ eta_jump_dp[pm]=exp(log_eta[i][j])-eta_cell[ip][jp];
+ dp_pos[pm](d)=grid_pos(d)
+ + (2*pm-1)*(intersection
+ (dist[i][j], dist_cell[ip][jp],
+ dist[i+(2*pm-1)*(d==0)][j+(2*pm-1)*(d==1)],h/2)
+ + h/2);
+ dp_pos[pm]((d+1)%2)=grid_pos((d+1)%2);
+ }
+ }
+ }
+
+ template <int ny_x, int ny_y>
+ void compute_xy(const int &d, const int &ii, const int &jj,
+ const double &log_eta_center,
+ const double log_etaY[][ny_y],
+ const double distX[][ny_x], const double distY[][ny_y],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
+ {
+ const int i(ii-(d==0)), j(jj-(d==1));
+
/* This has to be consistent with the checks for the second
derivatives. Otherwise, the stencil weights may get too large
for off-center points */
- intersect_sides[0]=(sign(dist[i][j])!=sign(dist[i+1][j])
- && sign(dist[i][j+1])!=sign(dist[i+1][j+1])
- && sign(dist[i][j])==sign(dist[i][j+1]));
- intersect_sides[1]=(sign(dist[i][j])!=sign(dist[i][j+1])
- && sign(dist[i+1][j])!=sign(dist[i+1][j+1])
- && sign(dist[i][j])==sign(dist[i+1][j]));
- int center(sign(center_dist));
+ int center(sign(distX[ii][jj]));
+
+
+ /* TODO: I do not quite like the explicit dependency on
+ intersect_dd. It seems like there could be cases where you do
+ want to use the side rather than corner interface locations,
+ but intersect_dd will not be true. */
+ intersect_sides[0][0]=(center!=sign(distY[i][j])
+ && center!=sign(distY[i][j+1])
+ && intersect_dd[0][0]);
+ intersect_sides[0][1]=(center!=sign(distY[i+1][j])
+ && center!=sign(distY[i+1][j+1])
+ && intersect_dd[0][1]);
+ intersect_sides[1][0]=(center!=sign(distY[i][j])
+ && center!=sign(distY[i+1][j])
+ && intersect_dd[1][0]);
+ intersect_sides[1][1]=(center!=sign(distY[i][j+1])
+ && center!=sign(distY[i+1][j+1])
+ && intersect_dd[1][1]);
+
+ /* TODO: fix for variable viscosity */
+
+ for(int dd=0;dd<2;++dd)
+ for(int pm=0;pm<2;++pm)
+ {
+ if(intersect_sides[dd][pm])
+ {
+ int ip=i+pm*(dd==0);
+ int jp=j+pm*(dd==1);
+
+ eta_jump_sides[dd][pm]=
+ exp(log_eta_center) - exp(log_etaY[ip][jp]);
+
+ if(d==dd)
+ {
+ sides_pos[dd][pm](dd)=grid_pos(dd) + (2*pm-1)
+ *(h/2 + intersection
+ (distX[ii][jj], dist_cell[ip][jp],
+ distX[ii+(2*pm-1)*(dd==0)][jj+(2*pm-1)*(dd==1)],
+ h/2));
+ sides_pos[dd][pm]((dd+1)%2)=grid_pos((dd+1)%2);
+ }
+ else
+ {
+ int ie=ii+pm*(dd==0);
+ int je=jj+pm*(dd==1);
+
+ sides_pos[dd][pm](dd)=grid_pos(dd) + (2*pm-1)
+ *(h/2 + intersection
+ (distX[ii][jj], dist_edge[ie][je],
+ distX[ii+(2*pm-1)*(dd==0)][jj+(2*pm-1)*(dd==1)],
+ h/2));
+ sides_pos[dd][pm]((dd+1)%2)=grid_pos((dd+1)%2);
+ }
+ }
+ }
+
+
for(int n0=0;n0<2;++n0)
for(int n1=0;n1<2;++n1)
{
- intersect_corner[n0][n1]=(center!=sign(dist[i+n0][j+n1]));
+ intersect_corner[n0][n1]=(center!=sign(distY[i+n0][j+n1]));
if(intersect_corner[n0][n1])
{
+ eta_jump_corner[n0][n1]=
+ exp(log_eta_center)-exp(log_etaY[i+n0][j+n1]);
+
FTensor::Tensor1<double,2> norm, tangent;
- norm(0)=dist[i+1][j]-dist[i][j];
- norm(1)=dist[i][j+1]-dist[i][j];
+ norm(0)=distY[i+1][j]-distY[i][j];
+ norm(1)=distY[i][j+1]-distY[i][j];
tangent(0)=-norm(1);
tangent(1)=norm(0);
@@ -65,9 +254,9 @@ public:
/* Top/Bottom */
corner_pos[n0][n1](0)=grid_pos(0)
+ sign(n0-0.5)*h*(0.5
- - std::fabs(dist[i+n0][j+n1]
- /(dist[i+(n0+1)%2][j+n1]
- - dist[i+n0][j+n1])));
+ - std::fabs(distY[i+n0][j+n1]
+ /(distY[i+(n0+1)%2][j+n1]
+ - distY[i+n0][j+n1])));
corner_pos[n0][n1](1)=grid_pos(1) + sign(n1-0.5)*h/2;
}
else
@@ -77,9 +266,9 @@ public:
corner_pos[n0][n1](1)=grid_pos(1)
+ sign(n1-0.5)*h*(0.5
- - std::fabs(dist[i+n0][j+n1]
- /(dist[i+n0][j+(n1+1)%2]
- - dist[i+n0][j+n1])));
+ - std::fabs(distY[i+n0][j+n1]
+ /(distY[i+n0][j+(n1+1)%2]
+ - distY[i+n0][j+n1])));
}
}
}
@@ -93,18 +282,21 @@ public:
&& intersect_corner[(n0+1)%2][(n1+1)%2];
if(intersect_anticorner[n0][n1])
{
+ eta_jump_anticorner[n0][n1]=
+ exp(log_eta_center)-exp(log_etaY[(n0+1)%2][n1]);
+
/* Only need to compute positions for the diagonal
points, since the center point has already been
calculated for interface_corner. */
- double delta=dist[i+n0][j+n1]-dist[i+(n0+1)%2][j+n1];
+ double delta=distY[i+n0][j+n1]-distY[i+(n0+1)%2][j+n1];
anticorner_pos[0](0)=grid_pos(0)
- +sign(n0-.5)*h*(0.5 - dist[i+n0][j+n1]/delta);
+ +sign(n0-.5)*h*(0.5 - distY[i+n0][j+n1]/delta);
anticorner_pos[0](1)=grid_pos(1)+sign(n1-.5)*h/2;
- delta=dist[i+n0][j+n1]-dist[i+n0][j+(n1+1)%2];
+ delta=distY[i+n0][j+n1]-distY[i+n0][j+(n1+1)%2];
anticorner_pos[1](0)=grid_pos(0)+sign(n0-.5)*h/2;
anticorner_pos[1](1)=grid_pos(1)
- +sign(n1-.5)*h*(0.5 - dist[i+n0][j+n1]/delta);
+ +sign(n1-.5)*h*(0.5 - distY[i+n0][j+n1]/delta);
std::cerr << "Anticorners not tested yet"
<< std::endl;
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/Interface/Interface.cxx
--- a/compute_coefficients/Interface/Interface.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/Interface/Interface.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -1,7 +1,13 @@
#include "../Interface.hxx"
Interface::Interface(const int &d, const int &i, const int &j,
- const double distx[Nx+1][Ny], const double disty[Nx][Ny+1])
+ const double log_etax[Nx+1][Ny],
+ const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
+ const double distx[Nx+1][Ny],
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
{
switch(d)
{
@@ -19,93 +25,16 @@ Interface::Interface(const int &d, const
if(d==0)
{
- intersect_dd[0]=sign(distx[i][j])!=sign(distx[i+1][j])
- || sign(distx[i][j])!=sign(distx[i-1][j]);
- // intersect_dd[0]=(i!=Nx-1 && sign(distx[i][j])!=sign(distx[i+1][j]))
- // || (i!=1 && sign(distx[i][j])!=sign(distx[i-1][j]));
- intersect_dd[1]=false;
- if(j!=Ny-1)
- {
- intersect_dd[1]=intersect_dd[1]
- || sign(distx[i][j])!=sign(distx[i][j+1]);
- }
- if(j!=0)
- {
- intersect_dd[1]=intersect_dd[1]
- || sign(distx[i][j-1])!=sign(distx[i][j]);
- }
-
- intersect_dp=(sign(disty[i][j] + disty[i][j+1])
- !=sign(disty[i-1][j] + disty[i-1][j+1]));
-
- /* TODO: This will have an error in the position of the interface
- of O(h^2). Need to correct for curvature. */
- double delta;
- if(intersect_dd[0])
- {
- dd_pos[0](1)=grid_pos(1);
- if(sign(distx[i][j])!=sign(distx[i+1][j]))
- delta=distx[i+1][j]-distx[i][j];
- else
- delta=distx[i][j]-distx[i-1][j];
-
- dd_pos[0](0)=grid_pos(0)-distx[i][j]*h/delta;
- }
- if(intersect_dd[1])
- {
- dd_pos[1](0)=grid_pos(0);
- if(j!=Ny-1 && sign(distx[i][j])!=sign(distx[i][j+1]))
- delta=distx[i][j+1]-distx[i][j];
- else
- delta=distx[i][j]-distx[i][j-1];
- dd_pos[1](1)=grid_pos(1)-distx[i][j]*h/delta;
- }
-
- compute_xy(i-1,j,distx[i][j],disty);
+ compute_dd(d,i,j,log_etax,distx,dist_cell,dist_edge);
+ compute_dp(d,i,j,log_etax,eta_cell,distx,dist_cell);
+ compute_xy(d,i,j,log_etax[i][j],log_etay,
+ distx,disty,dist_cell,dist_edge);
}
else
{
- intersect_dd[0]=false;
- if(i!=Nx-1)
- {
- intersect_dd[0]=intersect_dd[0]
- || sign(disty[i][j])!=sign(disty[i+1][j]);
- }
- if(i!=0)
- {
- intersect_dd[0]=intersect_dd[0]
- || sign(disty[i][j])!=sign(disty[i-1][j]);
- }
- intersect_dd[1]=sign(disty[i][j])!=sign(disty[i][j+1])
- || sign(disty[i][j-1])!=sign(disty[i][j]);
- // intersect_dd[1]=(j!=Ny-1 && sign(disty[i][j])!=sign(disty[i][j+1]))
- // || (j!=1 && sign(disty[i][j-1])!=sign(disty[i][j]));
-
- intersect_dp=(sign(distx[i][j] + distx[i+1][j])
- !=sign(distx[i][j-1] + distx[i+1][j-1]));
-
-
- double delta;
- if(intersect_dd[0])
- {
- dd_pos[0](1)=grid_pos(1);
- if(i!=Ny-1 && sign(disty[i][j])!=sign(disty[i+1][j]))
- delta=disty[i+1][j]-disty[i][j];
- else
- delta=disty[i][j]-disty[i-1][j];
-
- dd_pos[0](0)=grid_pos(0)-disty[i][j]*h/delta;
- }
- if(intersect_dd[1])
- {
- dd_pos[1](0)=grid_pos(0);
- if(sign(disty[i][j])!=sign(disty[i][j+1]))
- delta=disty[i][j+1]-disty[i][j];
- else
- delta=disty[i][j]-disty[i][j-1];
- dd_pos[1](1)=grid_pos(1)-disty[i][j]*h/delta;
- }
-
- compute_xy(i,j-1,disty[i][j],distx);
+ compute_dd(d,i,j,log_etay,disty,dist_cell,dist_edge);
+ compute_dp(d,i,j,log_etay,eta_cell,disty,dist_cell);
+ compute_xy(d,i,j,log_etay[i][j],log_etax,
+ disty,distx,dist_cell,dist_edge);
}
}
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/compute_Cp.cxx
--- a/compute_coefficients/compute_Cp.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/compute_Cp.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -3,6 +3,7 @@
#include "compute_v_on_interface.hxx"
#include "FTensor.hpp"
#include "sign.hxx"
+#include "Interface.hxx"
void compute_jumps(const double &eta_jump,
const FTensor::Tensor1<double,2> &v,
@@ -16,64 +17,79 @@ 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 eta_cell[Nx][Ny],
const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
const int &i, const int &j,
double &Cp)
{
Cp=0;
- bool intersects[]={sign(distx[i][j])!=sign(distx[i+1][j]),
- sign(disty[i][j])!=sign(disty[i][j+1])};
+ const int sgn(sign(dist_cell[i][j]));
+ bool intersects[][2]={{sign(distx[i][j])!=sgn,
+ sign(distx[i+1][j])!=sgn},
+ {sign(disty[i][j])!=sgn,
+ sign(disty[i][j+1])!=sgn}};
FTensor::Tensor1<double,2> pos((i+0.5)*h,(j+0.5)*h);
/* TODO fix for curved interfaces */
- FTensor::Tensor1<double,2> interface_pos[2];
- if(intersects[0])
+ FTensor::Tensor1<double,2> interface_pos[2][2];
+ for(int pm=0;pm<2;++pm)
{
- interface_pos[0](0)=i*h-distx[i][j]*h/(distx[i+1][j]-distx[i][j]);
- interface_pos[0](1)=(j+0.5)*h;
- }
- if(intersects[1])
- {
- interface_pos[1](0)=(i+0.5)*h;
- interface_pos[1](1)=j*h-disty[i][j]*h/(disty[i][j+1]-disty[i][j]);
+ if(intersects[0][pm])
+ {
+ interface_pos[0][pm](0)=pos(0)
+ + (1-2*pm)*Interface::intersection(distx[i+pm][j],dist_cell[i][j],
+ distx[i+1-pm][j],h/2);
+
+ interface_pos[0][pm](1)=pos(1);
+ }
+ if(intersects[1][pm])
+ {
+ interface_pos[1][pm](0)=pos(0);
+
+ interface_pos[1][pm](1)=pos(1)
+ + (1-2*pm)*Interface::intersection(disty[i][j+pm],dist_cell[i][j],
+ disty[i][j+1-pm],h/2);
+ }
}
for(int dd=0;dd<2;++dd)
- if(intersects[dd])
- {
- FTensor::Tensor1<double,2> v, dv, ddv;
- FTensor::Tensor1<int,2> dir(0,0);
- FTensor::Tensor2<double,2,2> nt_to_xy;
- const int sgn(sign(pos(dd)-interface_pos[dd](dd)));
- dir(dd)=sgn;
- compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface_pos[dd],dd,v,dv,ddv,nt_to_xy);
+ for(int pm=0;pm<2;++pm)
+ if(intersects[dd][pm])
+ {
+ FTensor::Tensor1<double,2> v, dv, ddv;
+ FTensor::Tensor1<int,2> dir(0,0);
+ FTensor::Tensor2<double,2,2> nt_to_xy;
+ const int sgn(sign(pos(dd)-interface_pos[dd][pm](dd)));
+ dir(dd)=sgn;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface_pos[dd][pm],dd,v,dv,ddv,nt_to_xy);
- FTensor::Tensor1<double,2> z_jump;
- FTensor::Tensor2<double,2,2> dz_jump;
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
- double eta_jump(sgn);
- switch(dd)
- {
- case 0:
- eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i][j]);
- break;
- case 1:
- eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j]);
- break;
- default:
- abort();
- break;
- }
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump);
+ double eta_jump;
+ switch(dd)
+ {
+ case 0:
+ eta_jump=eta_cell[i][j]-exp(log_etax[i+pm][j]);
+ break;
+ case 1:
+ eta_jump=eta_cell[i][j]-exp(log_etay[i][j+pm]);
+ break;
+ default:
+ abort();
+ break;
+ }
+ compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump);
- double dx=std::fabs(pos(dd)-interface_pos[dd](dd) - (h/2)*sgn);
+ double dx=std::fabs(pos(dd)-interface_pos[dd][pm](dd) - (h/2)*sgn);
- const FTensor::Index<'a',2> a;
- const FTensor::Index<'b',2> b;
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
- Cp-=dir(a)*(z_jump(a) - dx*dz_jump(a,b)*dir(b))/h;
- }
+ Cp-=dir(a)*(z_jump(a) - dx*dz_jump(a,b)*dir(b))/h;
+ }
}
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/compute_Cxyz.cxx
--- a/compute_coefficients/compute_Cxyz.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/compute_Cxyz.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -28,14 +28,18 @@ void compute_Cxyz(const double zx[Nx+1][
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1],
const int &d,
const int &i, const int &j,
double &C)
{
C=0;
- const Interface interface(d,i,j,distx,disty);
+ const Interface interface(d,i,j,log_etax,log_etay,eta_cell,
+ distx,disty,dist_cell,dist_edge);
if(!interface.intersects())
return;
@@ -48,116 +52,120 @@ void compute_Cxyz(const double zx[Nx+1][
FTensor::Tensor1<int,2> xyz(0,0);
xyz(d)=1;
bool try_corners(true);
- for(int dd=0;dd<2;++dd)
+ for(int pm=0;pm<2;++pm)
{
- if(interface.intersect_dd[dd])
+ for(int dd=0;dd<2;++dd)
{
- /* xyz= the component we are correcting, Cx, Cy, or Cz.
- dir=direction of 2nd derivative
- dir2, dir3=components of mixed derivative we are correcting.
- So for zy,xy with boundary in x, dir2=y, dir3=y. */
+ if(interface.intersect_dd[dd][pm])
+ {
+
+ /* xyz= the component we are correcting, Cx, Cy, or Cz.
+ dir=direction of 2nd derivative
+ dir2, dir3=components of mixed derivative we are correcting.
+ So for zy,xy with boundary in x, dir2=y, dir3=y. */
+ FTensor::Tensor1<double,2> v, dv, ddv;
+ FTensor::Tensor2<double,2,2> nt_to_xy;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface.dd_pos[dd][pm],d,
+ v,dv,ddv,nt_to_xy);
+
+ int sgn(sign(pos(dd) - interface.dd_pos[dd][pm](dd)));
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+ double p_jump;
+ FTensor::Tensor1<double,2> dp_jump;
+
+
+ compute_jumps(interface.eta_jump_dd[dd][pm],v,dv,ddv,nt_to_xy,
+ z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+
+ FTensor::Tensor1<double,2> dx(0,0);
+ dx(dd)=pos(dd) - interface.dd_pos[dd][pm](dd) - h*sgn;
+ double dz_dd_correction=
+ xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
+ + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
+
+ const int dz_dd_factor(d==dd ? 2 : 1);
+
+ C+=dz_dd_factor*dz_dd_correction;
+
+ // if(d==0 && dd==1 && i==9 && j==2)
+ // std::cout << __func__ << " "
+ // << d << " "
+ // // << pm << " "
+ // << i << " "
+ // << j << " "
+ // // << dx(0) << " "
+ // // << dx(1) << " "
+ // // << interface.dd_pos[dd][pm](0) << " "
+ // // << interface.dd_pos[dd][pm](1) << " "
+ // // << dv(0) << " "
+ // // << dv(1) << " "
+ // // << ddv(0) << " "
+ // // << ddv(1) << " "
+ // << interface.eta_jump_dd[dd][pm] << " "
+ // // << xyz(a)*(z_jump(a))/(h*h) << " "
+ // // << xyz(a)*(dz_jump(a,b)*dx(b))/(h*h) << " "
+ // << C << " "
+ // << "\n";
+
+ }
+ if(interface.intersect_sides[dd][pm])
+ {
+ FTensor::Tensor1<double,2> v, dv, ddv;
+ FTensor::Tensor2<double,2,2> nt_to_xy;
+ compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
+ interface.sides_pos[dd][pm],d,
+ v,dv,ddv,nt_to_xy);
+
+ int sgn(sign(pos(dd) - interface.sides_pos[dd][pm](dd)));
+ FTensor::Tensor1<double,2> z_jump;
+ FTensor::Tensor2<double,2,2> dz_jump;
+ FTensor::Tensor3_christof<double,2,2> ddz_jump;
+ double p_jump;
+ FTensor::Tensor1<double,2> dp_jump;
+
+ compute_jumps(interface.eta_jump_sides[dd][pm],v,dv,ddv,nt_to_xy,
+ z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
+
+ try_corners=false;
+ FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
+ dir2((d+1)%2)=1;
+ dir3((dd+1)%2)=1;
+
+ FTensor::Tensor1<double,2> dx(0,0);
+ dx(dd)=pos(dd) - interface.sides_pos[dd][pm](dd) - (h/2)*sgn;
+
+ C+=-sgn*(dz_jump(a,b)*dir2(a)*dir3(b)
+ + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
+ }
+ }
+ if(interface.intersect_dp[pm])
+ {
FTensor::Tensor1<double,2> v, dv, ddv;
FTensor::Tensor2<double,2,2> nt_to_xy;
compute_v_on_interface(zx,zy,log_etax,log_etay,distx,disty,
- interface.dd_pos[dd],d,
+ interface.dp_pos[pm],d,
v,dv,ddv,nt_to_xy);
- int sgn(sign(pos(dd) - interface.dd_pos[dd](dd)));
- FTensor::Tensor1<double,2> dx(0,0);
- dx(dd)=pos(dd) - interface.dd_pos[dd](dd) - h*sgn;
+ int sgn(sign(pos(d) - interface.dp_pos[pm](d)));
FTensor::Tensor1<double,2> z_jump;
FTensor::Tensor2<double,2,2> dz_jump;
FTensor::Tensor3_christof<double,2,2> ddz_jump;
double p_jump;
FTensor::Tensor1<double,2> dp_jump;
- double eta_jump(sgn);
- switch(d)
- {
- case 0:
- switch(dd)
- {
- case 0:
- eta_jump*=exp(log_etax[i+1][j])-exp(log_etax[i-1][j]);
- break;
- case 1:
- if(j==0)
- {
- eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j]);
- }
- else if(j==Ny-1)
- {
- eta_jump*=exp(log_etax[i][j])-exp(log_etax[i][j-1]);
- }
- else
- {
- eta_jump*=exp(log_etax[i][j+1])-exp(log_etax[i][j-1]);
- }
- break;
- default:
- abort();
- break;
- }
- break;
- case 1:
- switch(dd)
- {
- case 0:
- if(i==0)
- {
- eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i][j]);
- }
- else if(i==Nx-1)
- {
- eta_jump*=exp(log_etay[i][j])-exp(log_etay[i-1][j]);
- }
- else
- {
- eta_jump*=exp(log_etay[i+1][j])-exp(log_etay[i-1][j]);
- }
- break;
- case 1:
- eta_jump*=exp(log_etay[i][j+1])-exp(log_etay[i][j-1]);
- break;
- default:
- abort();
- break;
- }
- break;
- default:
- abort();
- break;
- }
+ compute_jumps(interface.eta_jump_dp[pm],v,dv,ddv,nt_to_xy,
+ z_jump,dz_jump,ddz_jump,p_jump,dp_jump);
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump,
- p_jump,dp_jump);
+ FTensor::Tensor1<double,2> dx(0,0);
+ dx(d)=pos(d) - interface.dp_pos[pm](d) - (h/2)*sgn;
- double dz_dd_correction=
- xyz(a)*(z_jump(a) + dz_jump(a,b)*dx(b)
- + ddz_jump(a,b,c)*dx(b)*dx(c)/2)/(h*h);
-
- const int dz_dd_factor(d==dd ? 2 : 1);
- C+=dz_dd_factor*dz_dd_correction;
-
- if(d==dd && interface.intersect_dp)
- {
- dx(dd)=pos(dd) - interface.dd_pos[dd](dd) - (h/2)*sgn;
- C+=sgn*(p_jump + dp_jump(a)*dx(a))/h;
- }
-
- if(interface.intersect_sides[dd])
- {
- try_corners=false;
- FTensor::Tensor1<int,2> dir2(0,0), dir3(0,0);
- dir2((d+1)%2)=1;
- dir3((dd+1)%2)=1;
- dx(dd)=pos(dd) - interface.dd_pos[dd](dd) - (h/2)*sgn;
-
- C+=-sgn*(dz_jump(a,b)*dir2(a)*dir3(b)
- + ddz_jump(a,b,c)*dir2(a)*dir3(b)*dx(c))/h;
- }
+ C+=sgn*(p_jump + dp_jump(a)*dx(a))/h;
}
}
+
for(int dd=0;dd<2;++dd)
for(int ee=0;ee<2;++ee)
if(interface.intersect_anticorner[dd][ee])
@@ -187,21 +195,8 @@ void compute_Cxyz(const double zx[Nx+1][
FTensor::Tensor2<double,2,2> dz_jump;
FTensor::Tensor3_christof<double,2,2> ddz_jump;
- double eta_jump;
- switch(d)
- {
- case 0:
- eta_jump=exp(log_etax[i][j])-exp(log_etay[i-1+dd][j+ee]);
- break;
-
- case 1:
- eta_jump=exp(log_etay[i][j])-exp(log_etax[i+dd][j-1+ee]);
- break;
- default:
- abort();
- break;
- }
- compute_jumps(eta_jump,v,dv,ddv,nt_to_xy,z_jump,dz_jump,ddz_jump);
+ compute_jumps(interface.eta_jump_corner[dd][ee],v,dv,ddv,nt_to_xy,
+ z_jump,dz_jump,ddz_jump);
FTensor::Tensor1<double,2> dx,yx_pos;
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/compute_Cxyz.hxx
--- a/compute_coefficients/compute_Cxyz.hxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/compute_Cxyz.hxx Fri Apr 13 15:46:54 2012 -0700
@@ -6,7 +6,10 @@ void compute_Cxyz(const double zx[Nx+1][
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double distx[Nx+1][Ny], const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1],
const int &d,
const int &i, const int &j,
double &C);
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/compute_coefficients.cxx
--- a/compute_coefficients/compute_coefficients.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/compute_coefficients.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -3,6 +3,8 @@
#include <vector>
#include <algorithm>
#include "FTensor.hpp"
+
+#include <iostream>
template< typename T, class Allocator >
void shrink_capacity(std::vector<T,Allocator>& v)
@@ -15,35 +17,46 @@ double simplified_Rx(const int &i, const
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double fx[Nx+1][Ny],
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1]);
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][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 eta_cell[Nx][Ny],
const double fx[Nx+1][Ny],
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1]);
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][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 eta_cell[Nx][Ny],
const double fx[Nx+1][Ny],
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1]);
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny]);
void compute_coefficients(const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][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],
@@ -76,14 +89,17 @@ void compute_coefficients(const double l
const int max_rp1(max_r+1);
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_rp1,0);di<=std::min(i+max_rp1,max_ij(0,0));++di)
- for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,max_ij(0,1));++dj)
+ for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(0,0));
+ ++di)
+ for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,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);
+ double R=simplified_Rx(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
+ distx,disty,dist_cell,dist_edge);
if(R!=0)
Cfx[0][i][j].push_back(Coefficient(di,dj,R));
@@ -91,39 +107,49 @@ void compute_coefficients(const double l
dRx[i][j]=R;
zx[di][dj]=0;
}
- for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(1,0));++di)
- for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,max_ij(1,1));++dj)
+ for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(1,0));
+ ++di)
+ for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,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);
+ double R=simplified_Rx(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
+ distx,disty,dist_cell,dist_edge);
if(R!=0)
Cfx[1][i][j].push_back(Coefficient(di,dj,R));
zy[di][dj]=0;
}
+
+ /* TODO: figure out whether I really need to do this and
+ whether this is a good place to do it? Maybe it should be
+ done after everything has been calculated? */
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_rp1,0);di<=std::min(i+max_rp1,max_ij(0,0));++di)
- for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,max_ij(0,1));++dj)
+ for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(0,0));
+ ++di)
+ for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,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);
+ double R=simplified_Ry(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
+ distx,disty,dist_cell,dist_edge);
if(R!=0)
Cfy[0][i][j].push_back(Coefficient(di,dj,R));
zx[di][dj]=0;
}
- for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(1,0));++di)
- for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,max_ij(1,1));++dj)
+ for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(1,0));
+ ++di)
+ for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,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);
+ double R=simplified_Ry(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
+ distx,disty,dist_cell,dist_edge);
if(R!=0)
Cfy[1][i][j].push_back(Coefficient(di,dj,R));
@@ -138,22 +164,26 @@ void compute_coefficients(const double l
for(int i=0;i<Nx;++i)
for(int j=0;j<Ny;++j)
{
- for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(0,0));++di)
- for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,max_ij(0,1));++dj)
+ for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(0,0));
+ ++di)
+ for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,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);
+ double R=simplified_Rp(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
+ distx,disty,dist_cell);
if(R!=0)
Cfp[0][i][j].push_back(Coefficient(di,dj,R));
zx[di][dj]=0;
}
- for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(1,0));++di)
- for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,max_ij(1,1));++dj)
+ for(int di=std::max(i-max_rp1,0);di<=std::min(i+max_rp1,max_ij(1,0));
+ ++di)
+ for(int dj=std::max(j-max_rp1,0);dj<=std::min(j+max_rp1,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);
+ double R=simplified_Rp(i,j,zx,zy,log_etax,log_etay,eta_cell,fx,fy,
+ distx,disty,dist_cell);
if(R!=0)
Cfp[1][i][j].push_back(Coefficient(di,dj,R));
zy[di][dj]=0;
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/simplified_Rp.cxx
--- a/compute_coefficients/simplified_Rp.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/simplified_Rp.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -6,8 +6,10 @@ extern void compute_Cp(const double zx[N
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
const int &i, const int &j,
double &Cp);
@@ -16,15 +18,17 @@ double simplified_Rp(const int &i, const
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double fx[Nx+1][Ny],
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1])
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny])
{
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);
+ compute_Cp(zx,zy,log_etax,log_etay,eta_cell,distx,disty,dist_cell,i,j,Cp);
return dzx_x + dzy_y + Cp;
}
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/simplified_Rx.cxx
--- a/compute_coefficients/simplified_Rx.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/simplified_Rx.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -8,10 +8,13 @@ double simplified_Rx(const int &i, const
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double fx[Nx+1][Ny],
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1])
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
{
double dzx_xx=
(zx[i-1][j] - 2*zx[i][j] + zx[i+1][j])/(h*h);
@@ -34,7 +37,8 @@ double simplified_Rx(const int &i, const
}
double Cx;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,0,i,j,Cx);
+ compute_Cxyz(zx,zy,log_etax,log_etay,eta_cell,distx,disty,dist_cell,dist_edge,
+ 0,i,j,Cx);
return 2*dzx_xx + dzx_yy + dzy_xy + Cx;
}
diff -r f5a306bd419e -r 8945333f6754 compute_coefficients/simplified_Ry.cxx
--- a/compute_coefficients/simplified_Ry.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/compute_coefficients/simplified_Ry.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -8,10 +8,13 @@ double simplified_Ry(const int &i, const
const double zy[Nx][Ny+1],
const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double fx[Nx+1][Ny],
const double fy[Nx][Ny+1],
const double distx[Nx+1][Ny],
- const double disty[Nx][Ny+1])
+ const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][Ny+1])
{
double dzy_yy=
(zy[i][j-1] - 2*zy[i][j] + zy[i][j+1])/(h*h);
@@ -34,7 +37,8 @@ double simplified_Ry(const int &i, const
}
double Cy;
- compute_Cxyz(zx,zy,log_etax,log_etay,distx,disty,1,i,j,Cy);
+ compute_Cxyz(zx,zy,log_etax,log_etay,eta_cell,distx,disty,dist_cell,dist_edge,
+ 1,i,j,Cy);
return 2*dzy_yy + dzy_xx + dzx_yx + Cy;
}
diff -r f5a306bd419e -r 8945333f6754 initial.cxx
--- a/initial.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/initial.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -9,8 +9,11 @@
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 eta_cell[Nx][Ny],
double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
- double distx[Nx+1][Ny], double disty[Nx][Ny+1])
+ double distx[Nx+1][Ny], double disty[Nx][Ny+1],
+ double dist_cell[Nx][Ny],
+ double dist_edge[Nx+1][Ny+1])
{
const double A(min_eta*(max_eta-min_eta)/(max_eta+min_eta));
std::complex<double> phi, psi, dphi, v;
@@ -52,9 +55,9 @@ void initial(const Model &model, double
gsl_permutation_free(perm);
- // const double theta(-0.5);
+ const double theta(-0.5);
// const double theta(pi);
- const double theta(0);
+ // 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);
@@ -64,6 +67,28 @@ void initial(const Model &model, double
for(int i=0;i<Nx+1;++i)
for(int j=0;j<Ny+1;++j)
{
+ {
+ FTensor::Tensor1<double,2> coord(i*h,j*h), xy;
+ xy(a)=rot(a,b)*(coord(b)-offset(b)) + offset(a);
+ // double x(xy(0)), y(xy(1)), r2(x*x+y*y);
+ double x(xy(0));
+
+ switch(model)
+ {
+ case Model::inclusion:
+ abort();
+ break;
+ case Model::solCx:
+ {
+ dist_edge[i][j]=x-middle;
+ }
+ break;
+ default:
+ abort();
+ break;
+ }
+ }
+
if(i<Nx && j<Ny)
{
FTensor::Tensor1<double,2> coord((i+0.5)*h,(j+0.5)*h), xy;
@@ -107,6 +132,8 @@ void initial(const Model &model, double
- sign*(kx*snh + csh - cs)*rho/(2*eta));
p[i][j]=(2*pi*eta*dvx_x - tau_xx)*std::cos(pi*y);
+ dist_cell[i][j]=x-middle;
+ eta_cell[i][j]=eta;
}
break;
default:
@@ -296,11 +323,13 @@ void initial(const Model &model, double
}
}
}
- write_vtk("p_initial",Nx,Ny,p);
- write_vtk("zx_initial",Nx+1,Ny,zx);
- write_vtk("zy_initial",Nx,Ny,zy);
- write_vtk("distx",Nx+1,Ny,distx);
- write_vtk("disty",Nx,Ny,disty);
- write_vtk("log_etax",Nx+1,Ny,log_etax);
- write_vtk("log_etay",Nx,Ny,log_etay);
+ write_vtk("p_initial",Nx,N,p);
+ write_vtk("zx_initial",Nx+1,N,zx);
+ write_vtk("zy_initial",Nx,N,zy);
+ write_vtk("distx",Nx+1,N,distx);
+ write_vtk("disty",Nx,N,disty);
+ write_vtk("dist_cell",Nx,N,dist_cell);
+ write_vtk("dist_edge",Nx+1,N,dist_edge);
+ write_vtk("log_etax",Nx+1,N,log_etax);
+ write_vtk("log_etay",Nx,N,log_etay);
}
diff -r f5a306bd419e -r 8945333f6754 main.cxx
--- a/main.cxx Sun Apr 08 17:37:14 2012 -0700
+++ b/main.cxx Fri Apr 13 15:46:54 2012 -0700
@@ -1,6 +1,7 @@
#include <iostream>
#include <fstream>
#include <sstream>
+#include <vector>
#include "Model.hxx"
#include "constants.hxx"
@@ -9,22 +10,19 @@
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 eta_cell[Nx][Ny],
double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1],
- double distx[Nx+1][Ny], double disty[Nx][Ny+1]);
-
-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 distx[Nx+1][Ny], double disty[Nx][Ny+1],
+ double dist_cell[Nx][Ny],
+ double dist_edge[Nx+1][Ny+1]);
void compute_coefficients(const double log_etax[Nx+1][Ny],
const double log_etay[Nx][Ny+1],
+ const double eta_cell[Nx][Ny],
const double distx[Nx+1][Ny],
const double disty[Nx][Ny+1],
+ const double dist_cell[Nx][Ny],
+ const double dist_edge[Nx+1][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],
@@ -40,8 +38,9 @@ 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], fx[Nx+1][Ny], fy[Nx][Ny+1],
- distx[Nx+1][Ny], disty[Nx][Ny+1], dRx[Nx+1][Ny], dRy[Nx][Ny+1];
+ eta_cell[Nx][Ny], p[Nx][Ny], dp[Nx][Ny], fx[Nx+1][Ny], fy[Nx][Ny+1],
+ distx[Nx+1][Ny], disty[Nx][Ny+1], dist_cell[Nx][Ny], dist_edge[Nx+1][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];
@@ -50,7 +49,7 @@ int main()
const bool jacobi(check_initial);
/* Initial conditions */
- const int n_sweeps(check_initial ? 1 : 40001);
+ const int n_sweeps(check_initial ? 1 : 1000001);
const double theta_mom=0.7;
const double theta_continuity=1.0;
const double tolerance=1.0e-6;
@@ -59,9 +58,11 @@ int main()
Model model(Model::solCx);
- initial(model,zx,zy,log_etax,log_etay,p,fx,fy,distx,disty);
+ initial(model,zx,zy,log_etax,log_etay,eta_cell,p,fx,fy,
+ distx,disty,dist_cell,dist_edge);
- compute_coefficients(log_etax,log_etay,distx,disty,Cfx,Cfy,Cfp,dRx,dRy);
+ compute_coefficients(log_etax,log_etay,eta_cell,distx,disty,dist_cell,
+ dist_edge,Cfx,Cfy,Cfp,dRx,dRy);
double Resid_p[Nx][Ny], Resid_x[Nx+1][Ny], Resid_y[Nx][Ny+1];
@@ -287,7 +288,8 @@ int main()
if(max_x_resid < tolerance
&& max_y_resid < tolerance
- && max_p_resid < tolerance)
+ )
+ // && max_p_resid < tolerance)
{
std::cout << "Solved "
<< sweep << " "
diff -r f5a306bd419e -r 8945333f6754 wscript
--- a/wscript Sun Apr 08 17:37:14 2012 -0700
+++ b/wscript Fri Apr 13 15:46:54 2012 -0700
@@ -22,6 +22,7 @@ def build(bld):
target = 'Gamr_disc',
# cxxflags = ['-std=c++0x','-g','-Wall','-Drestrict=','-DFTENSOR_DEBUG'],
cxxflags = ['-std=c++0x','-O3','-Wall','-Drestrict='],
+ # cxxflags = ['-std=c++0x','-O3','-Wall','-Drestrict=','-DNDEBUG'],
# cxxflags = ['-std=c++0x','-O3','-Wall','-Drestrict=','-pg'],
# linkflags = ['-pg'],
lib = ['boost_filesystem', 'boost_system', 'gsl', 'gslcblas'],
More information about the CIG-COMMITS
mailing list