[cig-commits] commit: First pass at computing interface velocity with matrix inversion. Only tested for x resid. Seems to work except at the boundaries.
Mercurial
hg at geodynamics.org
Mon Feb 13 20:08:05 PST 2012
changeset: 45:f2e7db30c967
user: Walter Landry <wlandry at caltech.edu>
date: Mon Feb 13 15:39:07 2012 -0800
files: compute_v_on_interface.cxx constants.hxx main.cxx
description:
First pass at computing interface velocity with matrix inversion. Only tested for x resid. Seems to work except at the boundaries.
diff -r 2ceeeafb95a9 -r f2e7db30c967 compute_v_on_interface.cxx
--- a/compute_v_on_interface.cxx Fri Feb 10 14:11:18 2012 -0800
+++ b/compute_v_on_interface.cxx Mon Feb 13 15:39:07 2012 -0800
@@ -1,162 +1,181 @@
#include <algorithm>
#include <cmath>
#include <iostream>
+#include "constants.hxx"
+#include <map>
+#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_spline.h>
-#include "constants.hxx"
+namespace {
+ template<int Ny>
+ void add_to_matrix(const int &n_terms, const int n, double z[][Ny],
+ double log_eta[][Ny], const int i, const int j,
+ const double &u0, const double &dxm, const double &dxp,
+ const double &dy,
+ double m[], double rhs[], const int sign=1)
+ {
+ rhs[n]+=sign*z[i][j]*exp(-log_eta[i][j]);
+ m[0+n*n_terms]+=sign*u0;
+ m[1+n*n_terms]+=sign*dy;
+ m[2+n*n_terms]+=sign*dy*dy/2;
+ m[3+n*n_terms]+=sign*dxm;
+ m[4+n*n_terms]+=sign*dxm*dy;
+ m[5+n*n_terms]+=sign*dxm*dxm/2;
+ m[6+n*n_terms]+=sign*dxp;
+ m[7+n*n_terms]+=sign*dxp*dy;
+ m[8+n*n_terms]+=sign*dxp*dxp;
+ }
+}
-template<class T>
-double interpolate_v(T z, T log_eta, const int &i, const int &j,
- const double &h, const double &x)
+template<int Ny>
+void interpolate_to_surface(double z[][Ny], double log_eta[][Ny],
+ const int &i, const int &j,
+ const double &delta_x, const bool &aligned,
+ double &v0, double &dv_y, double &dv_yy)
{
- double u0=z[i][j]*exp(-log_eta[i][j]);
- double up=z[i+1][j]*exp(-log_eta[i+1][j]);
- double um=z[i-1][j]*exp(-log_eta[i-1][j]);
+ const int n_terms(9);
+ /* Format of m is
- double du=(up-um)/(2*h);
- double ddu=(up - 2*u0 + um)/(h*h);
+ u0, u,y u,yy um,x um,xy um,xx up,x up,xy up,xx
- return u0 + x*du + x*x*ddu/2;
+ where y is parallel to the surface and x is
+ perpendicular */
+
+ double m[n_terms*n_terms];
+ double rhs[n_terms];
+
+ std::fill(m,m+n_terms*n_terms,0.0);
+ std::fill(rhs,rhs+n_terms,0.0);
+
+ double dx(delta_x);
+ int i_center(i), di(1);
+ if(dx>h/2)
+ {
+ di=-1;
+ i_center++;
+ dx=h-dx;
+ }
+
+ if(aligned)
+ {
+ /* Ordering is
+ i_center,j
+ i_center,j+1
+ i_center,j-1
+ i_center-1,j
+ i_center-2,j
+ i_center-1,j+1 - i_center+1,j-1
+
+ i_center+1,j
+ i_center+2,j
+ i_center+1,j+1 - i_center-1,j-1
+ */
+
+ add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,0,-di*dx,0,m,rhs);
+ add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,0,-di*dx,h,m,rhs);
+ add_to_matrix(n_terms,2,z,log_eta,i_center,j-1,1,0,-di*dx,-h,m,rhs);
+ add_to_matrix(n_terms,3,z,log_eta,i_center-di,j,1,0,-di*(dx+h),0,m,rhs);
+ add_to_matrix(n_terms,4,z,log_eta,i_center-2*di,j,1,0,-di*(dx+2*h),
+ 0,m,rhs);
+ add_to_matrix(n_terms,5,z,log_eta,i_center-di,j+1,1,0,-di*(dx+h),h,m,rhs);
+ add_to_matrix(n_terms,5,z,log_eta,i_center-di,j-1,1,0,-di*(dx+h),
+ -h,m,rhs,-1);
+
+ add_to_matrix(n_terms,6,z,log_eta,i_center+di,j,1,-di*(dx-h),0,0,m,rhs);
+ add_to_matrix(n_terms,7,z,log_eta,i_center+2*di,j,1,-di*(dx-2*h),0,
+ 0,m,rhs);
+ add_to_matrix(n_terms,8,z,log_eta,i_center+di,j+1,1,-di*(dx-h),0,h,m,rhs);
+ add_to_matrix(n_terms,8,z,log_eta,i_center+di,j-1,1,-di*(dx-h),0,
+ -h,m,rhs,-1);
+ }
+ else
+ {
+ /* Ordering is
+ i_center,j
+ i_center,j+1
+ i_center-1,j
+ i_center-1,j+1
+ i_center,j+2 + i_center,j-1
+ i_center-2,j + i_center-2,j+1
+
+ i_center+1,j
+ i_center+1,j+1
+ i_center+2,j + i_center+2,j+1
+ */
+
+ add_to_matrix(n_terms,0,z,log_eta,i_center,j,1,0,-di*dx,-h/2,m,rhs);
+ add_to_matrix(n_terms,1,z,log_eta,i_center,j+1,1,0,-di*dx,h/2,m,rhs);
+ add_to_matrix(n_terms,2,z,log_eta,i_center-di,j,1,0,-di*(dx+h),
+ -h/2,m,rhs);
+ add_to_matrix(n_terms,3,z,log_eta,i_center-di,j+1,1,0,-di*(dx+h),
+ h/2,m,rhs);
+ add_to_matrix(n_terms,4,z,log_eta,i_center,j+2,1,0,-di*dx,3*h/2,m,rhs);
+ add_to_matrix(n_terms,4,z,log_eta,i_center,j-1,1,0,-di*dx,-3*h/2,m,rhs);
+ add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j,1,0,-di*(dx+2*h),
+ -h/2,m,rhs);
+ add_to_matrix(n_terms,5,z,log_eta,i_center-2*di,j+1,1,0,-di*(dx+2*h),
+ h/2,m,rhs);
+
+ add_to_matrix(n_terms,6,z,log_eta,i_center+di,j,1,-di*(dx-h),0,-h/2,m,rhs);
+ add_to_matrix(n_terms,7,z,log_eta,i_center+di,j+1,1,-di*(dx-h),0,
+ h/2,m,rhs);
+ add_to_matrix(n_terms,8,z,log_eta,i_center+2*di,j,1,-di*(dx-2*h),0,
+ -h/2,m,rhs);
+ add_to_matrix(n_terms,8,z,log_eta,i_center+2*di,j+1,1,-di*(dx-2*h),0,
+ h/2,m,rhs);
+ }
+
+ // if(j==128)
+ // for(int a=0;a<9;++a)
+ // {
+ // std::cout << a << " "
+ // << rhs[a] << " ";
+ // for(int b=0;b<9;++b)
+ // std::cout << m[b+9*a] << " ";
+ // std::cout << "\n";
+ // }
+
+ gsl_matrix_view mv=gsl_matrix_view_array(m,n_terms,n_terms);
+ gsl_vector_view rhsv=gsl_vector_view_array(rhs,n_terms);
+ gsl_vector *result=gsl_vector_alloc(n_terms);
+ int s;
+ gsl_permutation *perm=gsl_permutation_alloc(n_terms);
+ gsl_linalg_LU_decomp(&mv.matrix,perm,&s);
+ gsl_linalg_LU_solve(&mv.matrix,perm,&rhsv.vector,result);
+
+ v0=gsl_vector_get(result,0);
+ dv_y=gsl_vector_get(result,1);
+ dv_yy=gsl_vector_get(result,2);
+ gsl_vector_free(result);
+ gsl_permutation_free(perm);
}
+
void compute_v_on_interface(double zx[Nx+1][Ny], double zy[Nx][Ny+1],
double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
- double vx[Ny+1], double dvx_y[Ny+1],
- double dvx_yy[Ny],
- double vy[Ny+1], double dvy_y[Ny],
- double dvy_yy[Ny+1])
+ const double &x, const int &j,
+ const int &xyz,
+ double &vx, double &vy,
+ double &dvx_y, double &dvy_y,
+ double &dvx_yy, double &dvy_yy)
{
- double y_x[Ny+1], y_y[Ny+1];
- for(int i=0;i<Ny+1;++i)
- y_x[i]=h*(i-0.5);
- for(int i=0;i<Ny+1;++i)
- y_y[i]=h*i;
+ int jx,jy;
+ if(xyz==0)
+ {
+ jx=jy=j;
+ }
+ else
+ {
+ jx=j-1;
+ jy=j;
+ }
- const double a1(h), a2(2*a1), b1(a1), b2(a2);
- double ux1_p[Ny], ux2_p[Ny], ux1_m[Ny], ux2_m[Ny],
- uy1_p[Ny+1], uy2_p[Ny+1], uy1_m[Ny+1], uy2_m[Ny+1];
+ /* vx */
+ int ix(x/h);
+ double dx_x(x-ix*h);
+ interpolate_to_surface(zx,log_etax,ix,jx,dx_x,xyz==0,vx,dvx_y,dvx_yy);
- {
- int i(middle/h);
- double dx=middle/h-i;
-
- /* vx, dvx_y, dvx_yy */
- for(int j=0;j<Ny;++j)
- {
- double low=interpolate_v(zx,log_etax,i-1,j,h,(1+dx)*h);
- double high=interpolate_v(zx,log_etax,i+2,j,h,(-2+dx)*h);
-
- vx[j+1]=low*(1-dx) + high*dx;
-
- ux1_p[j]=zx[i+1][j]*exp(-log_etax[i+1][j])*(1-dx)
- + zx[i+2][j]*exp(-log_etax[i+2][j])*dx;
- ux2_p[j]=zx[i+2][j]*exp(-log_etax[i+2][j])*(1-dx)
- + zx[i+3][j]*exp(-log_etax[i+3][j])*dx;
-
- ux1_m[j]=zx[i-1][j]*exp(-log_etax[i-1][j])*(1-dx)
- + zx[i][j]*exp(-log_etax[i][j])*dx;
- ux2_m[j]=zx[i-2][j]*exp(-log_etax[i-2][j])*(1-dx)
- + zx[i-1][j]*exp(-log_etax[i-1][j])*dx;
- }
- vx[0]=vx[Ny];
- }
- {
- int i(middle/h-0.5);
- double dx=middle/h-0.5-i;
-
- /* vy, dvy_y, dvy_yy */
- for(int j=0;j<Ny+1;++j)
- {
- double low=interpolate_v(zy,log_etay,i-1,j,h,h*(dx+1));
- double high=interpolate_v(zy,log_etay,i+2,j,h,h*(-2+dx));
-
- vy[j]=low*(1-dx) + high*dx;
-
- uy1_p[j]=zy[i+1][j]*exp(-log_etay[i+1][j])*(1-dx)
- + zy[i+2][j]*exp(-log_etay[i+2][j])*dx;
- uy2_p[j]=zy[i+2][j]*exp(-log_etay[i+2][j])*(1-dx)
- + zy[i+3][j]*exp(-log_etay[i+3][j])*dx;
-
- uy1_m[j]=zy[i-1][j]*exp(-log_etay[i-1][j])*(1-dx)
- + zy[i][j]*exp(-log_etay[i][j])*dx;
- uy2_m[j]=zy[i-2][j]*exp(-log_etay[i-2][j])*(1-dx)
- + zy[i-1][j]*exp(-log_etay[i-1][j])*dx;
- }
- }
-
- /* Iterate to find v on the boundary */
-
- const double relax=0.25;
- const double tolerance=1e-8;
- double error=0;
- do
- {
- error=0;
- /* Compute the derivatives from the current guess */
- /* dvx_dy, dvx_dyy */
- gsl_interp_accel *vx_accel=gsl_interp_accel_alloc ();
- gsl_spline *vx_spline=gsl_spline_alloc(gsl_interp_cspline_periodic, Ny+1);
- gsl_spline_init(vx_spline, y_x, vx, Ny+1);
-
- dvx_y[0]=0;
- dvx_y[Ny]=0;
- for(int i=1;i<Ny;++i)
- dvx_y[i]=gsl_spline_eval_deriv(vx_spline,y_y[i],vx_accel);
-
- for(int i=0;i<Ny;++i)
- dvx_yy[i]=gsl_spline_eval_deriv2(vx_spline,y_x[i+1],vx_accel);
-
- gsl_spline_free(vx_spline);
- gsl_interp_accel_free(vx_accel);
-
- /* dvy_dy, dvy_dyy */
- gsl_interp_accel *vy_accel=gsl_interp_accel_alloc ();
- gsl_spline *vy_spline=gsl_spline_alloc(gsl_interp_cspline_periodic, Ny+1);
- gsl_spline_init(vy_spline, y_y, vy, Ny+1);
-
- for(int i=0;i<Ny;++i)
- dvy_y[i]=gsl_spline_eval_deriv(vy_spline,y_x[i+1],vy_accel);
-
- for(int i=0;i<Ny;++i)
- dvy_yy[i]=gsl_spline_eval_deriv2(vy_spline,y_y[i],vy_accel);
-
- gsl_spline_free(vy_spline);
- gsl_interp_accel_free(vy_accel);
-
- /* Update the current guess based on the derivatives */
- for(int i=0;i<Ny;++i)
- {
- double zx_jump=-eta_jump*dvy_y[i];
- double temp=
- (min_eta*(a1*a2/(b2-b1))*(b2*b2*ux1_m[i] - b1*b1*ux2_m[i])
- + max_eta*(b1*b2/(a2-a1))*(a2*a2*ux1_p[i] - a1*a1*ux2_p[i]))
- /(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
- - a1*a2*b1*b2*zx_jump/(min_eta*a1*a2*(b1+b2)
- + max_eta*b1*b2*(a1+a2));
- if(vx[i+1]!=0 || temp!=0)
- error=
- std::max(error,std::fabs((vx[i+1]-temp)
- /(std::fabs(vx[i+1])+std::fabs(temp))));
- vx[i+1]=(1-relax)*vx[i+1] + relax*temp;
- }
- vx[0]=vx[Ny];
-
- for(int i=0;i<Ny+1;++i)
- {
- double zy_jump=-eta_jump*dvx_y[i];
- double temp=
- (min_eta*(a1*a2/(b2-b1))*(b2*b2*uy1_m[i] - b1*b1*uy2_m[i])
- + max_eta*(b1*b2/(a2-a1))*(a2*a2*uy1_p[i] - a1*a1*uy2_p[i]))
- /(min_eta*a1*a2*(b1+b2) + max_eta*b1*b2*(a1+a2))
- - a1*a2*b1*b2*zy_jump/(min_eta*a1*a2*(b1+b2)
- + max_eta*b1*b2*(a1+a2));
- if(vy[i]!=0 || temp!=0)
- error=
- std::max(error,std::fabs((vy[i]-temp)
- /(std::fabs(vy[i])+std::fabs(temp))));
-
- vy[i]=(1-relax)*vy[i] + relax*temp;
- }
- }
- while (error>tolerance);
+ /* vy */
+ int iy=x/h-0.5;
+ double dx_y=x-((iy+0.5)*h);
+ interpolate_to_surface(zy,log_etay,iy,jy,dx_y,xyz!=0,vy,dvy_y,dvy_yy);
}
diff -r 2ceeeafb95a9 -r f2e7db30c967 constants.hxx
--- a/constants.hxx Fri Feb 10 14:11:18 2012 -0800
+++ b/constants.hxx Mon Feb 13 15:39:07 2012 -0800
@@ -1,13 +1,14 @@ const int NN(256);
-const int NN(256);
+const int NN(512);
const int Nx(NN);
const int Ny(2*NN);
const double min_eta=1;
-const double max_eta=1e6;
+const double max_eta=1e4;
const double eta_jump=max_eta-min_eta;
#include <cmath>
const double log_max_eta=std::log(max_eta);
const double log_min_eta=std::log(min_eta);
-const double middle=25.00000001/64;
+const double middle=(25 + 0.00000001)/64;
+// const double middle=(25 + 32/Nx + 0.00000001)/64;
// const double middle=0.4;
const double h(1.0/Nx);
const double pi(atan(1.0)*4);
diff -r 2ceeeafb95a9 -r f2e7db30c967 main.cxx
--- a/main.cxx Fri Feb 10 14:11:18 2012 -0800
+++ b/main.cxx Mon Feb 13 15:39:07 2012 -0800
@@ -11,12 +11,12 @@ extern void initial(const Model &model,
double p[Nx][Ny], double fx[Nx+1][Ny], double fy[Nx][Ny+1]);
extern void compute_v_on_interface(double zx[Nx+1][Ny], double zy[Nx][Ny+1],
- double log_etax[Nx+1][Ny],
- double log_etay[Nx][Ny+1],
- double vx[Ny+1], double dvx_y[Ny+1],
- double dvx_yy[Ny],
- double vy[Ny+1], double dvy_y[Ny],
- double dvy_yy[Ny+1]);
+ double log_etax[Nx+1][Ny], double log_etay[Nx][Ny+1],
+ const double &x, const int &j,
+ const int &xyz,
+ double &vx, double &vy,
+ double &dvx_y, double &dvy_y,
+ double &dvx_yy, double &dvy_yy);
int main()
{
@@ -45,15 +45,6 @@ int main()
for(int iteration=0;iteration<max_iterations;++iteration)
{
- /* Compute the boundary jumps */
-
- double vx[Ny+1], dvx_y[Ny+1], dvx_yy[Ny],
- vy[Ny+1], dvy_y[Ny], dvy_yy[Ny+1];
-
- if(model==Model::solCx)
- compute_v_on_interface(zx,zy,log_etax,log_etay,vx,dvx_y,dvx_yy,
- vy,dvy_y,dvy_yy);
-
/* Smoothing sweeps */
// if(iteration==max_iterations-1)
@@ -145,19 +136,49 @@ int main()
if(x-h<middle && x+h>middle)
{
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ /* We waste some effort by computing
+ vx etc. once for each side,
+ though it does not change. */
+ compute_v_on_interface(zx,zy,log_etax,
+ log_etay,
+ middle,j,0,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
double dx=(x>middle) ? (middle-(x-h))
: (middle-(x+h));
- double zx_jump=eta_jump*vx[j+1];
- double dzx_x_jump=eta_jump*dvy_y[j];
- double p_jump=-2*dvy_y[j]*eta_jump;
- double dp_x_jump=2*dvx_yy[j]*eta_jump;
+ double zx_jump=eta_jump*vx;
+ double dzx_x_jump=eta_jump*dvy_y;
+ double p_jump=-2*dvy_y*eta_jump;
+ double dp_x_jump=2*dvx_yy*eta_jump;
double dzx_xx_jump=
- -dvx_yy[j]*eta_jump + dp_x_jump;
+ -dvx_yy*eta_jump + dp_x_jump;
double dzx_xx_correction=
-(zx_jump + dx*dzx_x_jump
+ dx*dx*dzx_xx_jump/2)/(h*h);
if(x>middle)
dzx_xx_correction*=-1;
+
+ // if(j==Ny/8)
+ // std::cout << "vx "
+ // << i << " "
+ // << j << " "
+ // << x << " "
+ // << middle << " "
+ // // << zx_jump << " "
+ // // << dzx_x_jump << " "
+ // // << dzx_xx_jump << " "
+ // // << dzx_xx << " "
+ // << (zx[i+3][j]-2*zx[i+2][j]+zx[i+1][j])/(h*h) << " "
+ // << (zx[i+2][j]-2*zx[i+1][j]+zx[i][j])/(h*h) << " "
+ // << (zx[i+1][j]-2*zx[i][j]+zx[i-1][j])/(h*h) << " "
+ // << (zx[i][j]-2*zx[i-1][j]+zx[i-2][j])/(h*h) << " "
+ // << (zx[i-1][j]-2*zx[i-2][j]+zx[i-3][j])/(h*h) << " "
+ // << dzx_xx+dzx_xx_correction << " "
+ // << "\n";
+
dzx_xx+=dzx_xx_correction;
@@ -165,10 +186,10 @@ int main()
{
dx=(x>middle) ? ((x-h/2)-middle)
: ((x+h/2)-middle);
+
dp_x-=(p_jump + dx*dp_x_jump)/h;
- dzy_xy-=
- eta_jump*(dvy_y[j] - dx*dvx_yy[j])/h;
+ dzy_xy-=eta_jump*(dvy_y - dx*dvx_yy)/h;
}
}
}
@@ -278,11 +299,22 @@ int main()
if(x-h<middle && x+h>middle)
{
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ /* We waste some effort by computing
+ vx etc. once for each side,
+ though it does not change. */
+ compute_v_on_interface(zx,zy,log_etax,
+ log_etay,
+ middle,j,1,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
double dx=(x>middle) ? (middle-(x-h))
: (middle-(x+h));
- double zy_jump=eta_jump*vy[j];
- double dzy_x_jump=-eta_jump*dvx_y[j];
- double dzy_xx_jump=-3*eta_jump*dvy_yy[j];
+ double zy_jump=eta_jump*vy;
+ double dzy_x_jump=-eta_jump*dvx_y;
+ double dzy_xx_jump=-3*eta_jump*dvy_yy;
double dzy_xx_correction=
-(zy_jump - dx*dzy_x_jump + dx*dx*dzy_xx_jump/2)/(h*h);
if(x>middle)
@@ -294,7 +326,7 @@ int main()
dx=(x>middle) ? ((x-h/2)-middle)
: ((x+h/2)-middle);
double dzx_yx_correction=
- eta_jump*(dvx_y[j] - dx*dvy_yy[j])/h;
+ eta_jump*(dvx_y - dx*dvy_yy)/h;
dzx_yx-=dzx_yx_correction;
// if(j==Ny/8)
// std::cout << "dzx_yx "
@@ -398,9 +430,19 @@ int main()
if(x+h/2>middle && x-h/2<middle)
{
+ double vx, dvx_y, dvx_yy, vy, dvy_y, dvy_yy;
+
+ /* We waste some effort by computing vx
+ etc. once for each side, though it does not
+ change. */
+ compute_v_on_interface(zx,zy,log_etax,log_etay,
+ middle,j,0,
+ vx,vy,dvx_y,dvy_y,
+ dvx_yy,dvy_yy);
+
double dx=(x>middle) ? (middle-(x-h/2)) : (middle-(x+h/2));
- double zx_jump=eta_jump*vx[j+1];
- double dzx_x_jump=eta_jump*dvy_y[j];
+ double zx_jump=eta_jump*vx;
+ double dzx_x_jump=eta_jump*dvy_y;
double dzx_x_correction=(zx_jump + dx*dzx_x_jump)/h;
dzx_x-=dzx_x_correction;
More information about the CIG-COMMITS
mailing list