NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
45-convdiff1d.cxx
Go to the documentation of this file.
1 
7 #include <iostream>
8 #include <limits>
9 #include <cassert>
10 #include <numcxx/numcxx.hxx>
11 #include <cmath>
12 #include <iostream>
13 
14 
15 #include "vtkfigFrame.h"
16 #include "vtkfigXYPlot.h"
17 
18 
19 
20 
21 inline double B(double x)
22 {
23  if (std::fabs(x)<1.0e-10) return 1.0;
24  return x/(std::exp(x)-1.0);
25 }
26 
27 void solve_expfit(double v, double D, numcxx::DArray1 &U)
28 {
29  int n=U.size();
30  auto h=1.0/(double)(n-1);
31  numcxx::DSparseMatrix A(n,n);
32  numcxx::DArray1 F(n);
33 
34  F=0;
35  U=0;
36  for (int k=0, l=1;k<n-1;k++,l++)
37  {
38  double g_kl=D* B(v*h/D);
39  double g_lk=D* B(-v*h/D);
40  A(k,k)+=g_kl/h;
41  A(k,l)-=g_kl/h;
42  A(l,l)+=g_lk/h;
43  A(l,k)-=g_lk/h;
44  }
45 
46  A(0,0)+=1.0e30;
47  A(n-1,n-1)+=1.0e30;
48  F(n-1)=1.0e30;
49  A.flush();
50 
51  numcxx::DSolverUMFPACK Solver(A);
52  Solver.update(A);
53  Solver.solve(U,F);
54 
55 }
56 
57 void solve_central(double v, double D, numcxx::DArray1 &U)
58 {
59  int n=U.size();
60  auto h=1.0/(double)(n-1);
61  numcxx::DSparseMatrix A(n,n);
62  numcxx::DArray1 F(n);
63 
64  F=0;
65  U=0;
66  double g_kl=D - 0.5*(v*h);
67  double g_lk=D + 0.5*(v*h);
68  for (int k=0, l=1;k<n-1;k++,l++)
69  {
70  A(k,k)+=g_kl/h;
71  A(k,l)-=g_kl/h;
72  A(l,l)+=g_lk/h;
73  A(l,k)-=g_lk/h;
74  }
75 
76  A(0,0)+=1.0e30;
77  A(n-1,n-1)+=1.0e30;
78  F(n-1)=1.0e30;
79  A.flush();
80 
81  numcxx::DSolverUMFPACK Solver(A);
82  Solver.update(A);
83  Solver.solve(U,F);
84 
85 }
86 
87 
88 void solve_upwind(double v, double D, numcxx::DArray1 &U)
89 {
90  int n=U.size();
91  auto h=1.0/(double)(n-1);
92  numcxx::DSparseMatrix A(n,n);
93  numcxx::DArray1 F(n);
94 
95  F=0;
96  U=0;
97  for (int k=0, l=1;k<n-1;k++,l++)
98  {
99  double g_kl=D;
100  double g_lk=D;
101  if (v<0) g_kl-=v*h;
102  else g_lk+=v*h;
103 
104  A(k,k)+=g_kl/h;
105  A(k,l)-=g_kl/h;
106  A(l,l)+=g_lk/h;
107  A(l,k)-=g_lk/h;
108  }
109 
110  A(0,0)+=1.0e30;
111  A(n-1,n-1)+=1.0e30;
112  F(n-1)=1.0e30;
113  A.flush();
114 
115  numcxx::DSolverUMFPACK Solver(A);
116  Solver.update(A);
117  Solver.solve(U,F);
118 
119 }
120 
121 
122 
123 
124 int main(void)
125 {
126 
127  int n=10;
128 
129  int maxref=10;
130  double D=0.01;
131  double v=1.0;
132 
133  auto frame=vtkfig::Frame::New();
134  frame->SetSize(800,800);
135  auto plot=vtkfig::XYPlot::New();
136  frame->AddFigure(plot);
137 
138  for (int iref=0;iref<maxref;iref++)
139  {
140  auto Uexp=numcxx::DArray1(n);
141  auto Uupw=numcxx::DArray1(n);
142  auto Ucnt=numcxx::DArray1(n);
143  auto X=numcxx::linspace(0.0,1.0,n);
144 
145 
146  solve_expfit(v,D,Uexp);
147  solve_central(v,D,Ucnt);
148  solve_upwind(v,D,Uupw);
149 
150 
151  frame->SetFrameTitle("n="+std::to_string(n));
152  plot->Clear();
153  plot->ShowGrid(false);
154 
155  plot->SetXRange(0,1);
156  plot->SetYRange(-1,1);
157 
158  plot->SetLegendPosition(0.3,0.7);
159  plot->SetPlotLineType("-");
160 
161 
162  plot->SetPlotColor(1,0,0);
163  plot->SetPlotLegend("exp");
164  plot->AddPlot(*X,Uexp);
165 
166  plot->SetPlotColor(0,1,0);
167  plot->SetPlotLegend("upw");
168  plot->AddPlot(*X,Uupw);
169 
170  plot->SetPlotColor(0,0,1);
171  plot->SetPlotLegend("cent");
172  plot->AddPlot(*X,Ucnt);
173 
174  frame->Interact();
175 
176  n*=2;
177  }
178 
179 }
180 
181 
Sparse matrix class using CRS storage scheme.
void flush()
Re-create the internal data structure in order to accomodated all newly created elements.
size_t size() const
Obtain size of array.
Definition: tarray.hxx:58
void solve(TArray< T > &Sol, const TArray< T > &Rhs)
Solve LU factorized system.
Main header of the library.
void solve_expfit(double v, double D, numcxx::DArray1 &U)
std::shared_ptr< TArray1< T > > linspace(T x0, T x1, int n)
Create array of $n$ equally spaced entries.
Definition: util.ixx:9
int main(void)
void solve_upwind(double v, double D, numcxx::DArray1 &U)
One dimensional array class.
Definition: tarray1.hxx:31
void solve_central(double v, double D, numcxx::DArray1 &U)
double B(double x)
Bridge class for using umfpack as solver for vmatrix.
TArray1< double > DArray1
Definition: numcxx.hxx:38
void update()
Perform actual computation of LU factorization.