NUMCXX  0.13.20181108
Numerical library for small projects and teaching purposes
Functions
48-convtest-fem-ilu.cxx File Reference
+ Include dependency graph for 48-convtest-fem-ilu.cxx:

Go to the source code of this file.

Functions

int main (void)
 

Function Documentation

int main ( void  )

Definition at line 25 of file 48-convtest-fem-ilu.cxx.

26 {
27  double maxref=9;
28  numcxx::Geometry Geometry;
29  Geometry.set_points({
30  {0,0},
31  {1,0},
32  {1,1},
33  {0,1},
34  {0.7,0.7},
35  });
36 
37  Geometry.set_bfaces({
38  {0,1},
39  {1,2},
40  {2,3},
41  {3,0}
42  });
43 
44 
45  Geometry.set_bfaceregions({1,2,3,4});
46 
47  Geometry.set_regionpoints({
48  {0.5,0.5}
49  });
50  Geometry.set_regionnumbers({1});
51 
52 
53  auto frame=vtkfig::Frame::New();
54  frame->SetSize(800,800);
55  frame->SetLayout(2,2);
56 
57  auto gridview=vtkfig::GridView::New();
58  frame->AddFigure(gridview,0);
59 
60  auto solview=vtkfig::ScalarView::New();
61  frame->AddFigure(solview,1);
62 
63  auto errplot=vtkfig::XYPlot::New();
64  frame->AddFigure(errplot,2);
65  errplot->SetXTitle("log10(h)");
66  errplot->SetYTitle("log10(error)");
67 
68  auto timeplot=vtkfig::XYPlot::New();
69  frame->AddFigure(timeplot,3);
70  timeplot->SetXTitle("log10(N)");
71  timeplot->SetYTitle("log10(time)");
72 
73 
74 
75 
76  std::vector<double> H;
77  std::vector<double> Cond;
78  std::vector<double> N;
79  std::vector<double> TAsm;
80  std::vector<double> TGen;
81  std::vector<double> TJac;
82  std::vector<double> TCG;
83 
84 
85  std::vector<double> L2Error;
86  std::vector<double> H1Error;
87 
88 
89  for (int iref=0;iref<maxref;iref++)
90  {
91  double t0;
92  double h_intended=0.5*pow(2.0,-iref);
93  double vol= h_intended*h_intended*0.25;
94 
95  Geometry.set_regionvolumes({vol});
96  t0=numcxx::cpu_clock();
97  numcxx::SimpleGrid grid(Geometry,"zpaAqDQ");
98  double t_gen=numcxx::cpu_clock()-t0;
99 
100  double hmin,hmax;
101  grid.calc_hminmax(hmin,hmax);
102  printf("h_intended: %8.3g h_max: %8.3g\n",h_intended,hmax);
103 
104 
105  numcxx::DArray1 bcfac(6);
106  numcxx::DArray1 bcval(6);
107  bcfac=0;
108  bcval=0;
109 
110  bcfac(1)=fem2d::Dirichlet;
111  bcfac(2)=fem2d::Dirichlet;
112  bcfac(3)=fem2d::Dirichlet;
113  bcfac(4)=fem2d::Dirichlet;
114  bcval(1)=0.0;
115  bcval(2)=0.0;
116  bcval(3)=0.0;
117  bcval(4)=0.0;
118 
119  auto nnodes=grid.npoints();
120  auto nodes=grid.get_points();
121 
122  numcxx::DArray1 source(nnodes);
123  numcxx::DArray1 exact(nnodes);
124  numcxx::DArray1 kappa(nnodes);
125  kappa=1;
126 
127  for (int i=0; i<nnodes;i++)
128  {
129  double x=nodes(i,0);
130  double y=nodes(i,1);
131  exact[i]=sin(M_PI*x)*sin(M_PI*y);
132  source[i]=2.0*M_PI*M_PI*exact[i];
133  }
134 
135  numcxx::DSparseMatrix SGlobal(nnodes,nnodes);
136  numcxx::DArray1 Rhs(nnodes);
137  numcxx::DArray1 Sol(nnodes);
138  numcxx::DArray1 Res(nnodes);
139  numcxx::DArray1 Upd(nnodes);
140 
141  t0=numcxx::cpu_clock();
142  fem2d::assemble_heat_problem(grid,bcfac,bcval,source,kappa,SGlobal, Rhs);
143  double t_asm=numcxx::cpu_clock()-t0;
144 
145  numcxx::DPreconILU Precon(SGlobal);
146 
147  t0=numcxx::cpu_clock();
148  Precon.update(SGlobal);
149 
150  t0=numcxx::cpu_clock();
151  Sol=0.0;
152  unsigned long max_iter=10000000;
153  double eps=1.0e-10;
154  double res0=numcxx::norm2(SGlobal*Sol-Rhs);
155  unsigned long nsteps=0;
156  for (;nsteps<max_iter;nsteps++)
157  {
158  Res=SGlobal*Sol-Rhs;
159  Precon.solve(Upd,Res);
160  Sol=Sol-Upd;
161  double norm=numcxx::norm2(Upd);
162  if (norm<1.0e-10) break;
163  }
164  double res=numcxx::norm2(SGlobal*Sol-Rhs);
165  double rho=pow(res/res0,1.0/nsteps);
166  double cond=(1.0+rho)/(1.0-rho);
167  printf("ILU contr: %g, per step: %g, kappa: %g\n", res/res0, rho,cond);
168  double t_jac=numcxx::cpu_clock()-t0;
169 
170 
171 
172  t0=numcxx::cpu_clock();
173  Sol=0.0;
174  res0=numcxx::norm2(SGlobal*Sol-Rhs);
175  int max_iter_cg=100000;
176  netlib::CG(SGlobal,Sol, Rhs,Precon,max_iter_cg,eps);
177  res=numcxx::norm2(SGlobal*Sol-Rhs);
178  rho=pow(res/res0,1.0/max_iter_cg);
179  printf("CG contr: %g, per step: %g\n", res/res0, rho);
180  double t_cg=numcxx::cpu_clock()-t0;
181 
182 
183 
184 
185 
186 
187  double l2error=fem2d::l2norm(grid,Sol-exact);
188  double h1error=fem2d::h1norm(grid,Sol-exact);
189 
190  printf("h: %8.3e l2: %8.3e h1: %8.3e\n",hmax,l2error,h1error);
191  H.push_back(log10(hmax));
192  H1Error.push_back(log10(h1error));
193  L2Error.push_back(log10(l2error));
194  Cond.push_back(log10(cond));
195 
196  printf("time/ms gen: %8.3f asm: %8.3f ilu: %8.3f cg: %8.3f\n",
197  t_gen*1000.0, t_asm*1000.0, t_jac*1000.0,t_cg*1000.0);
198 
199  N.push_back(log10(grid.npoints()));;
200  TAsm.push_back(log10(t_asm));
201  TJac.push_back(log10(t_jac));
202  TCG.push_back(log10(t_cg));
203  TGen.push_back(log10(t_gen));
204 
205 
206 #ifdef VTKFIG
207  auto griddata=numcxx::vtkfigDataSet(grid);
208  griddata->SetPointScalar(Sol ,"Sol");
209 
210  gridview->SetData(griddata);
211  solview->SetData(griddata,"Sol");
212 
213  frame->SetFrameTitle("iref="+std::to_string(iref));
214  if (iref>0)
215  {
216  double x0=-3;
217  double x1=0;
218  auto X=numcxx::DArray1{x0,x1};
219  auto OH1=numcxx::DArray1{x0,x1};
220  auto OH2=numcxx::DArray1{2*x0,2*x1};
221  auto OHM2=numcxx::DArray1{-2*x0-1,-2*x1-1};
222 
223  errplot->Clear();
224 
225 
226  errplot->SetLegendPosition(0.3,0.4);
227  errplot->SetXRange(x0,x1);
228  errplot->SetYRange(-5,5);
229 
230 
231  errplot->SetPlotLineType("-");
232  errplot->SetPlotMarkerType("o");
233  errplot->SetMarkerSize(2);
234 
235 
236  errplot->SetPlotColor(0.2,0.5,0.2);
237  errplot->SetPlotLegend("O(h**2)");
238  errplot->AddPlot(X,OH2);
239  errplot->SetPlotColor(0.5,0.2,0.2);
240  errplot->SetPlotLegend("O(h)");
241  errplot->AddPlot(X,OH1);
242 
243  errplot->SetPlotColor(0.0,0.5,0);
244  errplot->SetPlotLegend("L2");
245  errplot->AddPlot(H,L2Error);
246 
247  errplot->SetPlotLegend("H1");
248  errplot->SetPlotColor(0.5,0,0);
249  errplot->AddPlot(H,H1Error);
250 
251  errplot->SetPlotLegend("Cond");
252  errplot->SetPlotColor(0,0,0.5);
253  errplot->AddPlot(H,Cond);
254  errplot->SetPlotLegend("O(1/h**2)");
255  errplot->AddPlot(X,OHM2);
256 
257 
258  double n0=0;
259  double n1=7;
260  double n12=4;
261  double n11=6;
262  double t0=-5;
263  double t1=2;
264 
265  auto XN=numcxx::DArray1{n0,n1};
266  auto XN1=numcxx::DArray1{n0,n11};
267  auto XN2=numcxx::DArray1{n0,n12};
268  auto ON1=numcxx::DArray1{n0-6.0,n1-6.0};
269  auto ON2=numcxx::DArray1{2.0*n0-6, 2*n12-6};
270  auto ON32=numcxx::DArray1{1.5*n0-8.0, 1.5*n11-8.0};
271 
272 
273  timeplot->SetLegendPosition(0.3,0.6);
274  timeplot->SetLegendSize(0.15,0.3);
275  timeplot->SetXRange(n0,n1);
276  timeplot->SetYRange(t0,t1);
277 
278 
279  timeplot->Clear();
280  timeplot->SetPlotLineType("-");
281  timeplot->SetPlotMarkerType("o");
282  timeplot->SetMarkerSize(2);
283 
284 
285  // timeplot->SetPlotLegend("Gen");
286  // timeplot->SetPlotColor(0.5,0,0);
287  // timeplot->AddPlot(N,TGen);
288 
289  // timeplot->SetPlotLegend("Asm");
290  // timeplot->SetPlotColor(0.0,0.5,0);
291  // timeplot->AddPlot(N,TAsm);
292 
293  timeplot->SetPlotLegend("ILU");
294  timeplot->SetPlotColor(0.0,0.0,0.5);
295  timeplot->AddPlot(N,TJac);
296 
297  timeplot->SetPlotLegend("CG");
298  timeplot->SetPlotColor(0.0,0.5,0.5);
299  timeplot->AddPlot(N,TCG);
300 
301  timeplot->SetPlotLegend("O(N**2))");
302  timeplot->SetPlotColor(0.0,0.0,1.0);
303  timeplot->AddPlot(XN2,ON2);
304 
305  timeplot->SetPlotLegend("O(N**(3/2))");
306  timeplot->SetPlotColor(1.0,0.0,0.0);
307  timeplot->AddPlot(XN1,ON32);
308 
309  // timeplot->SetPlotLegend("O(N)");
310  // timeplot->SetPlotColor(1.0,0,0);
311  // timeplot->AddPlot(XN,ON1);;
312 
313 
314  }
315  frame->Interact();
316 #endif
317  }
318 
319 }
numcxx::Geometry
Class collecting data for the description of piecewise linear geometries.
Definition: geometry.hxx:17
numcxx::TSparseMatrix
Sparse matrix class using CRS storage scheme.
Definition: tsparsematrix.hxx:23
numcxx::Geometry::set_regionpoints
void set_regionpoints(const std::initializer_list< std::initializer_list< double >> &il)
Set member via intializer list.
Definition: geometry.hxx:32
numcxx::Geometry::set_regionvolumes
void set_regionvolumes(const std::initializer_list< double > &il)
Set member via intializer list.
Definition: geometry.hxx:38
netlib::CG
int CG(const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M, int &max_iter, Real &tol)
Iterative template routine – CG.
Definition: cg.hxx:28
numcxx::norm2
A::value_type norm2(const A &a)
Euklidean norm of array or expression.
Definition: util.ixx:47
fem2d::assemble_heat_problem
void assemble_heat_problem(const numcxx::SimpleGrid &Grid, const numcxx::DArray1 &BCfac, const numcxx::DArray1 &BCval, const numcxx::DArray1 &Source, const numcxx::DArray1 &Kappa, numcxx::DSparseMatrix &SGlobal, numcxx::DArray1 &Rhs)
Definition: fem2d.cxx:74
fem2d::h1norm
double h1norm(const numcxx::SimpleGrid &grid, const numcxx::DArray1 &u)
Definition: fem2d.cxx:302
numcxx::Geometry::set_points
void set_points(const std::initializer_list< std::initializer_list< double >> &il)
Set member via intializer list.
Definition: geometry.hxx:23
numcxx::Geometry::set_regionnumbers
void set_regionnumbers(const std::initializer_list< int > &il)
Set member via intializer list.
Definition: geometry.hxx:35
fem2d::Dirichlet
const double Dirichlet
BC value marking Dirichlet boundary condition.
Definition: fem2d.hxx:11
numcxx::TPreconILU
ILU preconditioner class.
Definition: tprecon-ilu.hxx:16
numcxx::Geometry::set_bfaces
void set_bfaces(const std::initializer_list< std::initializer_list< int >> &il)
Set member via intializer list.
Definition: geometry.hxx:26
numcxx::Geometry::set_bfaceregions
void set_bfaceregions(const std::initializer_list< int > &il)
Set member via intializer list.
Definition: geometry.hxx:29
numcxx::TArray1
One dimensional array class.
Definition: tarray1.hxx:31
numcxx::SimpleGrid
Class containing data for simple grid data structure.
Definition: simplegrid.hxx:19
numcxx::cpu_clock
double cpu_clock()
cpu time in seconds
Definition: util.ixx:101
numcxx::norm
A::value_type norm(const A &a)
Euklidean norm of array or expression.
Definition: util.hxx:54
fem2d::l2norm
double l2norm(const numcxx::SimpleGrid &grid, const numcxx::DArray1 &u)
Definition: fem2d.cxx:275

+ Here is the call graph for this function:

Generated on Thu Nov 8 2018 22:52:26 for NUMCXX by   doxygen 1.8.11