[cig-commits] commit: Compute v in a general manner.
Mercurial
hg at geodynamics.org
Wed Mar 21 16:14:20 PDT 2012
changeset: 105:f9b5c217e7a4
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Wed Mar 21 16:13:22 2012 -0700
files: compute_v_on_interface/compute_dv_dtt.cxx compute_v_on_interface/compute_v_on_interface.cxx compute_v_on_interface/compute_values.hxx
description:
Compute v in a general manner.
diff -r ce3129de86ff -r f9b5c217e7a4 compute_v_on_interface/compute_dv_dtt.cxx
--- a/compute_v_on_interface/compute_dv_dtt.cxx Wed Mar 21 04:26:49 2012 -0700
+++ b/compute_v_on_interface/compute_dv_dtt.cxx Wed Mar 21 16:13:22 2012 -0700
@@ -6,6 +6,7 @@
#include "vel.hxx"
#include "compute_2nd_derivs.hxx"
#include "compute_1st_derivs.hxx"
+#include "compute_values.hxx"
#include <tuple>
#include <algorithm>
#include <iostream>
@@ -42,6 +43,10 @@ void compute_dv_dtt(const double zx[Nx+1
dv_pm[0](a,b)=0;
dv_pm[1](a,b)=0;
+ FTensor::Tensor1<double,2> v_pm[2];
+ v_pm[0](a)=0;
+ v_pm[1](a)=0;
+
double eta_pm[2];
for(int d=0;d<2;++d)
@@ -55,6 +60,8 @@ void compute_dv_dtt(const double zx[Nx+1
d_valid(0).resize((2*max_r+1)*(2*max_r+1));
d_valid(1).resize((2*max_r+1)*(2*max_r+1));
+ std::vector<Valid> valid((2*max_r+1)*(2*max_r+1));
+
int max_x(Nx-d), max_y(Ny+d-1);
double dx(d*h/2), dy((1-d)*h/2);
@@ -65,18 +72,20 @@ void compute_dv_dtt(const double zx[Nx+1
for(int j=std::max(startj-max_r,0);j<=std::min(startj+max_r,max_y);++j)
{
FTensor::Tensor1<double,2> p(i*h+dx,j*h+dy), p_d(i*h+dx+h/2,j*h+dy),
- dd_diff, d_diff;
+ dd_diff, d_diff, diff;
dd_diff(a)=p(a)-pos(a);
d_diff(a)=p_d(a)-pos(a);
if(d==0)
{
+ /* ddv */
if(i>0 && i<max_x)
dd_valid(0,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])
&& sign(distx[i][j])==sign(distx[i-1][j]),
i,j,dd_diff);
+ /* dv */
if(i<max_x)
{
d_diff(a)-=length*norm(a)*sign(distx[i][j]);
@@ -85,9 +94,14 @@ void compute_dv_dtt(const double zx[Nx+1
sign(distx[i][j])==sign(distx[i+1][j]),
i,j,d_diff);
}
+ /* v */
+ diff(a)=dd_diff(a)-1.5*length*norm(a)*sign(distx[i][j]);
+ valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(distx[i][j]),true,i,j,diff);
}
else
{
+ /* ddv */
bool v_dd;
if(i==0)
v_dd=(sign(disty[i][j])==sign(disty[i+1][j]));
@@ -100,6 +114,7 @@ void compute_dv_dtt(const double zx[Nx+1
dd_valid(0,0)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(disty[i][j]),v_dd,i,j,dd_diff);
+ /* dv */
bool v_d;
/* This bothers me a little. What if the interface
lies between the point and the boundary? Shouldn't
@@ -112,6 +127,11 @@ void compute_dv_dtt(const double zx[Nx+1
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);
+
+ /* v */
+ diff(a)=dd_diff(a)-1.5*length*norm(a)*sign(disty[i][j]);
+ valid[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
+ Valid(sign(disty[i][j]),true,i,j,diff);
}
}
@@ -126,6 +146,7 @@ void compute_dv_dtt(const double zx[Nx+1
if(d==0)
{
+ /* ddv */
bool v;
if(j==0)
v=sign(distx[i][j])==sign(distx[i][j+1]);
@@ -137,6 +158,7 @@ void compute_dv_dtt(const double zx[Nx+1
dd_valid(1,1)[i-starti+max_r + (2*max_r+1)*(j-startj+max_r)]=
Valid(sign(distx[i][j]),v,i,j,dd_diff);
+ /* dv */
bool v_d;
if(j==max_y)
v_d=true;
@@ -149,6 +171,7 @@ void compute_dv_dtt(const double zx[Nx+1
}
else
{
+ /* ddv */
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]),
@@ -156,6 +179,7 @@ void compute_dv_dtt(const double zx[Nx+1
&& sign(disty[i][j])==sign(disty[i][j-1]),
i,j,dd_diff);
+ /* dv */
if(j<max_y)
{
d_diff(a)-=length*norm(a)*sign(disty[i][j]);
@@ -207,19 +231,72 @@ void compute_dv_dtt(const double zx[Nx+1
}
}
+ /* Sort valid vectors */
+
std::sort(dd_valid(0,0).begin(),dd_valid(0,0).end());
std::sort(dd_valid(0,1).begin(),dd_valid(0,1).end());
std::sort(dd_valid(1,1).begin(),dd_valid(1,1).end());
std::sort(d_valid(0).begin(),d_valid(0).end());
std::sort(d_valid(1).begin(),d_valid(1).end());
+ std::sort(valid.begin(),valid.end());
FTensor::Tensor3_christof<bool,2,2> is_set(false,false,false,false,
false,false);
+ /* Compute everything */
+ /* v */
+
+ std::vector<Valid> values[2];
+ for(auto &V: valid)
+ {
+ const int pm=(V.sign==-1 ? 0 : 1);
+ switch(values[pm].size())
+ {
+ case 0:
+ case 1:
+ case 2:
+ values[pm].push_back(V);
+ break;
+ case 6:
+ break;
+ default:
+ int x(0),y(0);
+ std::vector<int> x_coords={V.i};
+ std::vector<int> y_coords={V.j};
+ for(auto &value: values[pm])
+ {
+ if(V.i==value.i)
+ ++x;
+ if(V.j==value.j)
+ ++y;
+ if(values[pm].size()==5)
+ {
+ if(find(x_coords.begin(),x_coords.end(),value.i)==x_coords.end())
+ x_coords.push_back(value.i);
+ if(find(y_coords.begin(),y_coords.end(),value.j)==y_coords.end())
+ y_coords.push_back(value.j);
+ }
+ }
+ if(x<3 && y<3 && !(values[pm].size()==5 &&
+ (x_coords.size()<3 || y_coords.size()<3)))
+ values[pm].push_back(V);
+ if(values[pm].size()==6)
+ {
+ if(d==0)
+ compute_values(0,values[pm],Nx+1,zx,log_etax,
+ norm,pm==0 ? -length : length,v_pm[pm]);
+ else
+ compute_values(1,values[pm],Nx,zy,log_etay,
+ norm,pm==0 ? -length : length,v_pm[pm]);
+ }
+ break;
+ }
+ }
+
for(int d0=0;d0<2;++d0)
{
+ /* First derivative dv */
std::vector<Valid> derivs[2];
- /* First derivative */
for(auto &V: d_valid(d0))
{
if(V.valid)
@@ -251,7 +328,7 @@ void compute_dv_dtt(const double zx[Nx+1
}
}
}
- /* Second derivative */
+ /* Second derivative ddv */
for(int d1=d0;d1<2;++d1)
{
for(auto &V: dd_valid(d0,d1))
@@ -303,5 +380,18 @@ void compute_dv_dtt(const double zx[Nx+1
dv(n)=dv_xy(a)*norm(a);
dv(t)=dv_xy(a)*tangent(a);
+ /* v */
+ FTensor::Tensor1<double,2> v_xy(0,0);
+ for(int d=0;d<2;++d)
+ v_xy(a)+=v_pm[d](a)*eta_pm[d];
+
+ FTensor::Tensor1<double,2> dz_n_jump;
+ dz_n_jump(a)=-eta_jump*dv(b)*ddel(a,b);
+
+ v_xy(a)-=2*h*dz_n_jump(a)/3;
+ v_xy(a)/=(eta_pm[0] + eta_pm[1]);
+
+ v(n)=v_xy(a)*norm(a);
+ v(t)=v_xy(a)*tangent(a);
}
diff -r ce3129de86ff -r f9b5c217e7a4 compute_v_on_interface/compute_v_on_interface.cxx
--- a/compute_v_on_interface/compute_v_on_interface.cxx Wed Mar 21 04:26:49 2012 -0700
+++ b/compute_v_on_interface/compute_v_on_interface.cxx Wed Mar 21 16:13:22 2012 -0700
@@ -6,45 +6,6 @@
#include "FTensor.hpp"
#include "../compute_v_on_interface.hxx"
#include "vel.hxx"
-
-template<int ny>
-double compute_v(const double z[][ny],
- const double log_eta[][ny],
- const double &jump,
- const double &dx,
- const int i, const int j)
-{
- /* v=v0 + dv*x + ddv*h*h/2, where x is centered on v_p2. The
- distance between points is 1 here, which simplifies the
- expressions. */
-
- double v_p1(vel(z,log_eta,i+1,j)), v_p2(vel(z,log_eta,i+2,j)),
- v_p3(vel(z,log_eta,i+3,j));
- double vp0(v_p2), dvp((v_p3-v_p1)/2), ddvp(v_p3 - 2*v_p2 + v_p1);
- double v_p(vp0 - (1-dx)*dvp + (1-dx)*(1-dx)*ddvp/2),
- v_pp(vp0 + dx*dvp + dx*dx*ddvp/2);
-
- double v_m0(vel(z,log_eta,i,j)), v_m1(vel(z,log_eta,i-1,j)),
- v_m2(vel(z,log_eta,i-2,j));
- double vm0(v_m1), dvm((v_m0-v_m2)/2), ddvm(v_m2 - 2*v_m1 + v_m0);
- double v_m(vm0 + dx*dvm + dx*dx*ddvm/2),
- v_mm(vm0 - (1-dx)*dvm + (1-dx)*(1-dx)*ddvm/2);
-
- return (4*(max_eta*v_p + min_eta*v_m) - (max_eta*v_pp + min_eta*v_mm)
- - 2*h*jump)
- /(3*(max_eta+min_eta));
-}
-
-/* For now, hard code the solCx answer */
-double dzx_x_jump(const double &dvy_y)
-{
- return -eta_jump*dvy_y;
-}
-
-double dzy_x_jump(const double &dvx_y)
-{
- return -eta_jump*dvx_y;
-}
void compute_dv_dtt(const double zx[Nx+1][Ny],
const double zy[Nx][Ny+1],
@@ -121,9 +82,6 @@ void compute_v_on_interface(const double
compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
i,j,ddv,ddv_dv_temp,dv,dv_dv_temp,v,v_dv_temp);
- v(0)=compute_v(zx,log_etax,dzx_x_jump(dv(1)),dx_i,ix,j);
- v(1)=0;
-
double eta=exp(log_etax[i][j]);
double x_i(i*h);
double delta_x=(x_i>middle) ? (middle-(x_i-h))
@@ -165,10 +123,7 @@ void compute_v_on_interface(const double
FTensor::Tensor1<double,2> ddv_dv_temp,dv_dv_temp,v_temp,v_dv_temp;
compute_dv_dtt(zx,zy,log_etax,log_etay,distx,disty,pos,norm,tangent,
- i,j,ddv,ddv_dv_temp,dv,dv_dv_temp,v_temp,v_dv_temp);
-
- v(1)=compute_v(zy,log_etay,dzy_x_jump(dv(0)),dx_j,iy,j);
- v(0)=0;
+ i,j,ddv,ddv_dv_temp,dv,dv_dv_temp,v,v_dv_temp);
double eta=exp(log_etay[i][j]);
double x_i((i+0.5)*h);
diff -r ce3129de86ff -r f9b5c217e7a4 compute_v_on_interface/compute_values.hxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_v_on_interface/compute_values.hxx Wed Mar 21 16:13:22 2012 -0700
@@ -0,0 +1,58 @@
+#ifndef GAMR_COMPUTE_VALUES_HXX
+#define GAMR_COMPUTE_VALUES_HXX
+
+#include <cassert>
+#include "Valid.hxx"
+#include <gsl/gsl_linalg.h>
+
+template<int ny>
+void compute_values(const int &d, const std::vector<Valid> &values,
+ const int &nx,
+ const double z[][ny],
+ const double log_eta[][ny],
+ const FTensor::Tensor1<double,2> &norm,
+ const double &signed_length,
+ FTensor::Tensor1<double,2> &v_pm)
+{
+ /* Create and invert the Vandermonde matrix */
+ assert(values.size()==6);
+ double m[6*6], rhs[6];
+ for(int i=0;i<6;++i)
+ {
+ rhs[i]=vel(z,log_eta,values[i].i,values[i].j);
+
+ m[0+6*i]=1;
+ m[1+6*i]=values[i].diff(0);
+ m[2+6*i]=values[i].diff(1);
+ m[3+6*i]=values[i].diff(0)*values[i].diff(0)/2;
+ m[4+6*i]=values[i].diff(1)*values[i].diff(1)/2;
+ m[5+6*i]=values[i].diff(0)*values[i].diff(1);
+ }
+
+ gsl_matrix_view mv=gsl_matrix_view_array(m,6,6);
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,6);
+ gsl_vector *x=gsl_vector_alloc(6);
+ int s;
+ gsl_permutation *perm=gsl_permutation_alloc(6);
+ gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
+ gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,x);
+
+ double v=gsl_vector_get(x,0);
+
+ FTensor::Tensor1<double,2> dv_xy;
+ const FTensor::Index<'a',2> a;
+ const FTensor::Index<'b',2> b;
+ dv_xy(0)=gsl_vector_get(x,1);
+ dv_xy(1)=gsl_vector_get(x,2);
+ double dv=dv_xy(a)*norm(a)*signed_length/2;
+
+ FTensor::Tensor2_symmetric<double,2> ddv_xy;
+ ddv_xy(0,0)=gsl_vector_get(x,3);
+ ddv_xy(1,1)=gsl_vector_get(x,4);
+ ddv_xy(0,1)=gsl_vector_get(x,5);
+ double ddv=ddv_xy(a,b)*norm(a)*norm(b)*signed_length*signed_length/4;
+
+ v_pm(d)=(3*v - 5*dv + 1.5*ddv)/3;
+}
+
+#endif
More information about the CIG-COMMITS
mailing list