45-convdiff1d.cxx File Reference

Include dependency graph for 45-convdiff1d.cxx:

Go to the source code of this file.

Functions
double	B (double x)
void	solve_expfit (double v, double D, numcxx::DArray1 &U)
void	solve_central (double v, double D, numcxx::DArray1 &U)
void	solve_upwind (double v, double D, numcxx::DArray1 &U)
int	main (void)

Function Documentation

double B ( double  x )

inline

Examples:

10-numcxx-expressions.cxx, and 45-convdiff1d.cxx.

Definition at line 21 of file 45-convdiff1d.cxx.

{
  if (std::fabs(x)<1.0e-10) return 1.0;
  return x/(std::exp(x)-1.0);
}

void solve_expfit	(	double	v,
		double	D,
		numcxx::DArray1 &	U
	)

Examples:

45-convdiff1d.cxx.

Definition at line 27 of file 45-convdiff1d.cxx.

{
  int n=U.size();
  auto h=1.0/(double)(n-1);
  numcxx::DSparseMatrix A(n,n);
  numcxx::DArray1 F(n);
    
  F=0;
  U=0;
  for (int k=0, l=1;k<n-1;k++,l++)
  {
    double g_kl=D* B(v*h/D);
    double g_lk=D* B(-v*h/D);
    A(k,k)+=g_kl/h;
    A(k,l)-=g_kl/h;
    A(l,l)+=g_lk/h;
    A(l,k)-=g_lk/h;
  }

  A(0,0)+=1.0e30;
  A(n-1,n-1)+=1.0e30;
  F(n-1)=1.0e30;
  A.flush();

  numcxx::DSolverUMFPACK Solver(A);
  Solver.update(A);
  Solver.solve(U,F);

}

Here is the call graph for this function:

void solve_central	(	double	v,
		double	D,
		numcxx::DArray1 &	U
	)

Examples:

45-convdiff1d.cxx.

Definition at line 57 of file 45-convdiff1d.cxx.

{
  int n=U.size();
  auto h=1.0/(double)(n-1);
  numcxx::DSparseMatrix A(n,n);
  numcxx::DArray1 F(n);
    
  F=0;
  U=0;
  double g_kl=D - 0.5*(v*h);
  double g_lk=D + 0.5*(v*h);
  for (int k=0, l=1;k<n-1;k++,l++)
  {
    A(k,k)+=g_kl/h;
    A(k,l)-=g_kl/h;
    A(l,l)+=g_lk/h;
    A(l,k)-=g_lk/h;
  }
  
  A(0,0)+=1.0e30;
  A(n-1,n-1)+=1.0e30;
  F(n-1)=1.0e30;
  A.flush();
  
  numcxx::DSolverUMFPACK Solver(A);
  Solver.update(A);
  Solver.solve(U,F);

}

Here is the call graph for this function:

void solve_upwind	(	double	v,
		double	D,
		numcxx::DArray1 &	U
	)

Examples:

45-convdiff1d.cxx.

Definition at line 88 of file 45-convdiff1d.cxx.

{
  int n=U.size();
  auto h=1.0/(double)(n-1);
  numcxx::DSparseMatrix A(n,n);
  numcxx::DArray1 F(n);
    
  F=0;
  U=0;
  for (int k=0, l=1;k<n-1;k++,l++)
  {
    double g_kl=D;
    double g_lk=D;
    if (v<0) g_kl-=v*h;
    else  g_lk+=v*h;

    A(k,k)+=g_kl/h;
    A(k,l)-=g_kl/h;
    A(l,l)+=g_lk/h;
    A(l,k)-=g_lk/h;
  }

  A(0,0)+=1.0e30;
  A(n-1,n-1)+=1.0e30;
  F(n-1)=1.0e30;
  A.flush();

  numcxx::DSolverUMFPACK Solver(A);
  Solver.update(A);
  Solver.solve(U,F);

}

Here is the call graph for this function:

int main

(

void

)

Examples:

45-convdiff1d.cxx.

Definition at line 124 of file 45-convdiff1d.cxx.

{

  int n=10;
  
  int maxref=10;
  double D=0.01;
  double v=1.0;
  
  auto frame=vtkfig::Frame::New();
  frame->SetSize(800,800);
  auto plot=vtkfig::XYPlot::New();
  frame->AddFigure(plot);
  
  for (int iref=0;iref<maxref;iref++)
  {
    auto Uexp=numcxx::DArray1(n);
    auto Uupw=numcxx::DArray1(n);
    auto Ucnt=numcxx::DArray1(n);
    auto X=numcxx::linspace(0.0,1.0,n);


    solve_expfit(v,D,Uexp);
    solve_central(v,D,Ucnt);
    solve_upwind(v,D,Uupw);


    frame->SetFrameTitle("n="+std::to_string(n));
    plot->Clear();
    plot->ShowGrid(false);

    plot->SetXRange(0,1);
    plot->SetYRange(-1,1);

    plot->SetLegendPosition(0.3,0.7);
    plot->SetPlotLineType("-");


    plot->SetPlotColor(1,0,0);
    plot->SetPlotLegend("exp");
    plot->AddPlot(*X,Uexp);

    plot->SetPlotColor(0,1,0);
    plot->SetPlotLegend("upw");
    plot->AddPlot(*X,Uupw);

    plot->SetPlotColor(0,0,1);
    plot->SetPlotLegend("cent");
    plot->AddPlot(*X,Ucnt);

    frame->Interact();

    n*=2;
  }

}

Here is the call graph for this function: