[cig-commits] commit: Compute dv_dt using general method.
Mercurial
hg at geodynamics.org
Wed Mar 21 03:47:56 PDT 2012
changeset: 103:1ca7a2a82be9
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Wed Mar 21 03:47:38 2012 -0700
files: compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx
description:
Compute dv_dt using general method.
diff -r 018b068cbb9f -r 1ca7a2a82be9 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Tue Mar 20 16:31:16 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Wed Mar 21 03:47:38 2012 -0700
@@ -21,9 +21,13 @@ void compute_dv_dtt(const double zx[Nx+1
const FTensor::Tensor1<double,2> &tangent,
const int &ii, const int &jj,
FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor1<double,2> &ddv_dv)
+ FTensor::Tensor1<double,2> &ddv_dv,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &dv_dv)
{
const int max_r=4;
+
+ double length=h/std::max(std::abs(norm(0)),std::abs(norm(1)));
const FTensor::Index<'a',2> a;
const FTensor::Index<'b',2> b;
@@ -63,7 +67,6 @@ void compute_dv_dtt(const double zx[Nx+1
dd_diff(a)=p(a)-pos(a);
d_diff(a)=p_d(a)-pos(a);
-
if(d==0)
{
if(i>0 && i<max_x)
@@ -73,11 +76,13 @@ void compute_dv_dtt(const double zx[Nx+1
&& sign(distx[i][j])==sign(distx[i-1][j]),
i,j,dd_diff);
if(i<max_x)
- d_valid(0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(distx[i][j]),
- sign(distx[i][j])==sign(distx[i+1][j]),
- i,j,d_diff);
-
+ {
+ d_diff(a)-=length*norm(a)*sign(distx[i][j]);
+ d_valid(0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),
+ sign(distx[i][j])==sign(distx[i+1][j]),
+ i,j,d_diff);
+ }
}
else
{
@@ -102,6 +107,7 @@ void compute_dv_dtt(const double zx[Nx+1
else
v_d=(sign(disty[i][j])==sign(disty[i+1][j]));
+ d_diff(a)-=length*norm(a)*sign(disty[i][j]);
d_valid(0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(disty[i][j]),v_d,i,j,d_diff);
}
@@ -134,6 +140,8 @@ void compute_dv_dtt(const double zx[Nx+1
v_d=true;
else
v_d=sign(distx[i][j])==sign(distx[i][j+1]);
+
+ d_diff(a)-=length*norm(a)*sign(distx[i][j]);
d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(distx[i][j]),v_d,i,j,d_diff);
}
@@ -141,13 +149,19 @@ void compute_dv_dtt(const double zx[Nx+1
{
if(j>0 && j<max_y)
dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),sign(disty[i][j])==sign(disty[i][j+1])
+ Valid(sign(disty[i][j]),
+ sign(disty[i][j])==sign(disty[i][j+1])
&& sign(disty[i][j])==sign(disty[i][j-1]),
i,j,dd_diff);
+
if(j<max_y)
- d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
- Valid(sign(disty[i][j]),sign(disty[i][j])==sign(disty[i][j+1]),
+ {
+ d_diff(a)-=length*norm(a)*sign(disty[i][j]);
+ d_valid(1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),
+ sign(disty[i][j])==sign(disty[i][j+1]),
i,j,d_diff);
+ }
}
}
@@ -257,16 +271,35 @@ void compute_dv_dtt(const double zx[Nx+1
}
}
}
+
+ /* ddv */
double temp(0);
FTensor::Tensor1<double,2> ddv_xy(0,0);
+ FTensor::Number<0> n;
+ FTensor::Number<1> t;
for(int d=0;d<2;++d)
{
temp+=eta_pm[d]*eta_pm[d];
ddv_xy(a)+=ddv_pm[d](a,b,c)*tangent(b)*tangent(c)*eta_pm[d]*eta_pm[d];
}
ddv_xy(a)/=temp*h*h;
+ ddv(n)=ddv_xy(a)*norm(a);
+ ddv(t)=ddv_xy(a)*tangent(a);
- ddv(0)=ddv_xy(a)*norm(a);
- ddv(1)=ddv_xy(a)*tangent(a);
+ /* dv */
+ FTensor::Tensor1<double,2> dv_xy(0,0);
+ for(int d=0;d<2;++d)
+ dv_xy(a)+=dv_pm[d](a,b)*tangent(b)*eta_pm[d];
+
+ FTensor::Tensor2_symmetric<int,2> ddel(0,1,0);
+ FTensor::Tensor1<double,2> ddz_nt_jump;
+ double eta_jump(eta_pm[1]-eta_pm[0]);
+ ddz_nt_jump(a)=-eta_jump*ddv(b)*ddel(a,b);
+
+ dv_xy(a)-=h*h*ddz_nt_jump(a);
+ dv_xy(a)/=(eta_pm[0]+eta_pm[1])*h;
+ dv(n)=dv_xy(a)*norm(a);
+ dv(t)=dv_xy(a)*tangent(a);
+
}
diff -r 018b068cbb9f -r 1ca7a2a82be9 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Tue Mar 20 16:31:16 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Wed Mar 21 03:47:38 2012 -0700
@@ -6,25 +6,6 @@
#include "FTensor.hpp"
#include "../compute_v_on_interface.hxx"
#include "vel.hxx"
-
-template<int ny>
-double compute_dv_y(const double z[][ny],
- const double log_eta[][ny],
- const double &jump,
- const double &dx,
- const int i, const int j)
-{
- double dv_y_p1=(vel(z,log_eta,i+1,j+1) - vel(z,log_eta,i+1,j))/h;
- double dv_y_p2=(vel(z,log_eta,i+2,j+1) - vel(z,log_eta,i+2,j))/h;
- double dv_y_p=(1-dx)*dv_y_p1 + dx*dv_y_p2;
-
- double dv_y_m0=(vel(z,log_eta,i,j+1) - vel(z,log_eta,i,j))/h;
- double dv_y_m1=(vel(z,log_eta,i-1,j+1) - vel(z,log_eta,i-1,j))/h;
- double dv_y_m=(1-dx)*dv_y_m1 + dx*dv_y_m0;
-
- return (max_eta*dv_y_p + min_eta*dv_y_m - h*jump)
- /(max_eta + min_eta);
-}
template<int ny>
double compute_v(const double z[][ny],
@@ -86,7 +67,9 @@ void compute_dv_dtt(const double zx[Nx+1
const FTensor::Tensor1<double,2> &tangent,
const int &i, const int &j,
FTensor::Tensor1<double,2> &ddv,
- FTensor::Tensor1<double,2> &ddv_dv);
+ FTensor::Tensor1<double,2> &ddv_dv,
+ FTensor::Tensor1<double,2> &dv,
+ FTensor::Tensor1<double,2> &dv_dv);
void compute_v_on_interface(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -142,23 +125,12 @@ void compute_v_on_interface(const double
tangent(0)=-norm(1);
tangent(1)=norm(0);
- FTensor::Tensor1<double,2> ddv_dv_temp;
+ FTensor::Tensor1<double,2> ddv_dv_temp,dv_dv_temp;
compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,
- norm,tangent,i,j,ddv,ddv_dv_temp);
+ norm,tangent,i,j,ddv,ddv_dv_temp,dv,dv_dv_temp);
- // double length=1/std::max(std::abs(norm(0)),std::abs(norm(1)));
-
- // FTensor::Tensor1<double,2> pos_p, pos_pp, pos_m, pos_mm;
- // pos_p(a)=pos(a)+length*norm(a);
- // pos_pp(a)=pos(a)+2*norm(a)*length;
- // pos_m(a)=pos(a)-length*norm(a);
- // pos_mm(a)=pos(a)-2*norm(a)*length;
-
- dv(1)=compute_dv_y(zy,log_etay,dzy_yx_jump(ddv(0)),dx_j,iy,j);
v(0)=compute_v(zx,log_etax,dzx_x_jump(dv(1)),dx_i,ix,j);
-
v(1)=0;
- dv(0)=0;
double eta=exp(log_etax[i][j]);
double x_i(i*h);
@@ -199,16 +171,12 @@ void compute_v_on_interface(const double
tangent(0)=-norm(1);
tangent(1)=norm(0);
- FTensor::Tensor1<double,2> ddv_temp, ddv_dv_temp;
+ FTensor::Tensor1<double,2> ddv_temp, ddv_dv_temp,dv_dv_temp;
compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,
- norm,tangent,i,j,ddv,ddv_dv_temp);
+ norm,tangent,i,j,ddv,ddv_dv_temp,dv,dv_dv_temp);
- dv(0)=compute_dv_y(zx,log_etax,dzx_yx_jump(ddv(1)),dx_i,ix,j-1);
v(1)=compute_v(zy,log_etay,dzy_x_jump(dv(0)),dx_j,iy,j);
-
v(0)=0;
- dv(1)=0;
- ddv(0)=0;
double eta=exp(log_etay[i][j]);
double x_i((i+0.5)*h);
More information about the CIG-COMMITS
mailing list