OperatorsTest.cpp 12.5 KB
Newer Older
1
#include <amdis/AMDiS.hpp>
2
#include <amdis/ContextGeometry.hpp>
3
#include <amdis/ProblemStat.hpp>
4
#include <amdis/LocalOperators.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
5
#include <amdis/algorithm/ForEach.hpp>
6
7
#include <amdis/localoperators/ConvectionDiffusionOperator.hpp>
#include <amdis/localoperators/StokesOperator.hpp>
8
9
10

using namespace AMDiS;

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
using T = double;

template <class LOp>
using IsOperatorList = decltype(std::declval<LOp>().empty());

// Test the functionality of a local operator
template <class LOp, class LocalView, class MatVec, class... Nodes>
void test_local_operator(LOp&& lop, LocalView const& localView, MatVec& matVec,
                         Nodes const&... nodes)
{
  auto const& gridView = localView.globalBasis().gridView();
  auto const& element = localView.element();
  auto geometry = element.geometry();

  [[maybe_unused]] auto context = [&]{
    if constexpr(Dune::Std::is_detected_v<IsOperatorList,LOp>)
Praetorius, Simon's avatar
Praetorius, Simon committed
27
      return GlobalContext<TYPEOF(gridView)>{gridView, element, geometry};
28
    else
Praetorius, Simon's avatar
Praetorius, Simon committed
29
      return ContextGeometry<TYPEOF(element)>{element, element, geometry};
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
  }();

  lop.bind(element);
  if constexpr(!Dune::Std::is_detected_v<IsOperatorList,LOp>)
    lop.assemble(context, nodes..., matVec);
  lop.unbind();
}


// Test the functionality of an operator
template <class Op, class LocalView, class... Args>
void test_operator(Op&& op, LocalView const& localView, Args&&... args)
{
  op.update(localView.globalBasis().gridView());

  auto lop = localOperator(op);
  test_local_operator(lop, localView, FWD(args)...);
}


// Test a hierarchic list of operators
template <class OpList, class LocalView, class... Args>
void test_operator_list(OpList& opList, LocalView const& localView, Args&&... args)
{
  Recursive::forEach(opList, [&](auto& op) {
    test_operator(op, localView, FWD(args)...);
  });

58
  auto lopList = localOperators(opList);
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
  Recursive::forEach(lopList, [&](auto& lop) {
    test_local_operator(lop, localView, FWD(args)...);
  });
}


// Test a tree-vector of operators
template <class OpList, class LocalView, class ElementVector>
void test_operator_list_vec(OpList& opList, LocalView const& localView,
                            ElementVector& vec)
{
  Traversal::forEachNode(localView.treeCache(), [&](auto const& node, auto&& tp) {
    test_operator_list(opList[tp].element_, localView, vec, node);
  });

74
  auto lopList = localOperators(opList);
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
  Traversal::forEachNode(localView.treeCache(), [&](auto const& node, auto&& tp) {
    Recursive::forEach(lopList[tp].element_, [&](auto& lop) {
      test_local_operator(lop, localView, vec, node);
    });
  });
}


// Test a tree-matrix of operators
template <class OpList, class LocalView, class ElementMatrix>
void test_operator_list_mat(OpList& opList, LocalView const& localView,
                            ElementMatrix& mat)
{
  Traversal::forEachNode(localView.treeCache(), [&](auto const& rowNode, auto&& r) {
    Traversal::forEachNode(localView.treeCache(), [&](auto const& colNode, auto&& c) {
      test_operator_list(opList[r][c].element_, localView, mat, rowNode, colNode);
    });
  });

94
  auto lopList = localOperators(opList);
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
  Traversal::forEachNode(localView.treeCache(), [&](auto const& rowNode, auto&& r) {
    Traversal::forEachNode(localView.treeCache(), [&](auto const& colNode, auto&& c) {
      Recursive::forEach(lopList[r][c].element_, [&](auto& lop) {
        test_local_operator(lop, localView, mat, rowNode, colNode);
      });
    });
  });
}


// Test the copy and move constructed operators
template <class Op, class Test>
void test_construction(Op const& op, Test test)
{
  // copy construction
  Op op2(op);
  test(op2);

  auto op3 = op2;
  test(op3);

  // move construction
  Op op4(std::move(op3));
  test(op4);

  auto op5 = std::move(op4);
  test(op5);
}


// Test a vector operator
template <class Op, class Basis, class Path>
void test(Op& op, Basis const& basis, Path path)
{
  using GridView = typename Basis::GridView;
  GridView gridView = basis.gridView();

  auto localView = basis.localView();
  localView.bind(*gridView.template begin<0>());

  auto const& tree = localView.treeCache();
  auto const& node = Dune::TypeTree::child(tree, path);

  using ElementVector = Dune::DynamicVector<T>;
  ElementVector vec(localView.size());

  using Element = typename GridView::template Codim<0>::Entity;
  auto a_op = makeOperator<Element>(op, gridView);
  test_operator(a_op, localView, vec, node);
  test_construction(a_op, [&](auto& op)
  {
    test_operator(op, localView, vec, node);
  });

  using Traits = Impl::OperatorTraits<GridView,Element,ElementVector>;

  // Construct a type-erasure Operator from a concrete operator
  using Node = TYPEOF(node);
  using OpBase = Operator<Traits,Node>;
  OpBase b_op(a_op);
  test_operator(b_op, localView, vec, node);
  test_construction(b_op, [&](auto& op)
  {
    test_operator(op, localView, vec, node);
  });

  VectorOperators<Basis,ElementVector> opList;
  opList[path].push(tag::element_operator<Element>{}, a_op);
  test_operator_list(opList, localView, vec, node);
  test_operator_list_vec(opList, localView, vec);

  OpBase c_op(std::move(a_op));
  test_operator(c_op, localView, vec, node);
}


// Test a matrix operator
template <class Op, class Basis, class Row, class Col>
void test(Op& op, Basis const& basis, Row row, Col col)
{
  using GridView = typename Basis::GridView;
  GridView gridView = basis.gridView();

  auto localView = basis.localView();
  localView.bind(*gridView.template begin<0>());

  auto const& tree = localView.treeCache();
  auto const& rowNode = Dune::TypeTree::child(tree, row);
  auto const& colNode = Dune::TypeTree::child(tree, col);

  using ElementMatrix = Dune::DynamicMatrix<T>;
  ElementMatrix mat(localView.size(), localView.size());

  using Element = typename GridView::template Codim<0>::Entity;
  auto a_op = makeOperator<Element>(op, gridView);
  test_operator(a_op, localView, mat, rowNode, colNode);
  test_construction(a_op, [&](auto& op)
  {
    test_operator(op, localView, mat, rowNode, colNode);
  });

  using Traits = Impl::OperatorTraits<GridView,Element,ElementMatrix>;

  // Construct a type-erasure Operator from a concrete operator
  using RowNode = TYPEOF(rowNode);
  using ColNode = TYPEOF(colNode);
  using OpBase = Operator<Traits,RowNode,ColNode>;
  OpBase b_op(a_op);
  test_operator(b_op, localView, mat, rowNode, colNode);
  test_construction(b_op, [&](auto& op)
  {
    test_operator(op, localView, mat, rowNode, colNode);
  });

  MatrixOperators<Basis,Basis,ElementMatrix> opList;
  opList[row][col].push(tag::element_operator<Element>{}, a_op);
  test_operator_list(opList, localView, mat, rowNode, colNode);
  test_operator_list_mat(opList, localView, mat);
}

template <class Problem>
void test_problem_operators(Problem& prob)
{
  auto localView = prob.globalBasis()->localView();
  localView.bind(*prob.gridView().template begin<0>());

  FlatVector<T> elementVec(localView.size());
  FlatMatrix<T> elementMat(localView.size(), localView.size());

  { // test problem operators
    test_operator_list_vec(prob.rhsVector()->operators(), localView, elementVec);
    test_operator_list_mat(prob.systemMatrix()->operators(), localView, elementMat);

    test_construction(prob.rhsVector()->operators(), [&](auto& op) {
      test_operator_list_vec(op, localView, elementVec);
    });
    test_construction(prob.systemMatrix()->operators(), [&](auto& op) {
      test_operator_list_mat(op, localView, elementMat);
    });
  }
}

237
238
239

int main(int argc, char** argv)
{
240
  Environment env(argc, argv);
241
242
243

  using namespace Dune::Indices;

244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
  // use T as coordinate type
  using HostGrid = Dune::YaspGrid<2, Dune::EquidistantCoordinates<T,2>>;
  HostGrid hostGrid({T(1), T(1)}, {8,8});
  AdaptiveGrid<HostGrid> grid(hostGrid);

  // use T as range type for basis
  using namespace Dune::Functions::BasisFactory;
#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
  GlobalBasis basis{grid.leafGridView(), composite(power<2>(lagrange<2>(), flatLexicographic()), lagrange<1>(), flatLexicographic())};
#else
  GlobalBasis basis{grid.leafGridView(), composite(power<2>(lagrange<2,T>(), flatLexicographic()), lagrange<1,T>(), flatLexicographic())};
#endif
  using Basis = decltype(basis);

  // use T as coefficient type
  using Param = DefaultProblemTraits<Basis, T>;


  ProblemStat<Param> prob("ellipt", grid, basis);
263
264
  prob.initialize(INIT_ALL);

265
266
  auto _v = makeTreePath(_0);
  auto _p = makeTreePath(_1);
267

268
269
270
271
272
273
274
275
276
277
  FieldVector<T,1> scalar = 1.0;
  FieldVector<T,2> vec; vec = 1.0;
  FieldMatrix<T,2,2> mat; mat = 1.0;

  auto op0a = sot(scalar);
  auto op0b = zot(scalar);
  test(op0a, *prob.globalBasis(), _p, _p);
  prob.addMatrixOperator(op0a, _p, _p);
  test(op0b, *prob.globalBasis(), _p);
  prob.addVectorOperator(op0b, _p);
278
279

  auto op1 = makeOperator(tag::divtestvec{}, scalar);
280
  test(op1, *prob.globalBasis(), _v);
281
282
283
  prob.addVectorOperator(op1, _v);

  auto op2 = makeOperator(tag::divtestvec_trial{}, scalar);
284
  test(op2, *prob.globalBasis(), _v, _p);
285
286
287
  prob.addMatrixOperator(op2, _v, _p);

  auto op3 = makeOperator(tag::gradtest{}, vec);
288
  test(op3, *prob.globalBasis(), _p);
289
290
291
  prob.addVectorOperator(op3, _p);

  auto op4 = makeOperator(tag::gradtest_trial{}, vec);
292
  test(op4, *prob.globalBasis(), _p, _p);
293
294
295
296
  prob.addMatrixOperator(op4, _p, _p);

  auto op5a = makeOperator(tag::gradtest_trialvec{}, scalar);
  // auto op5b = makeOperator(tag::gradtest_trialvec{}, mat); // TODO: needs to be implemented
297
  test(op5a, *prob.globalBasis(), _p, _v);
298
299
300
  prob.addMatrixOperator(op5a, _p, _v);

  auto op6 = makeOperator(tag::partialtest{0}, scalar);
301
  test(op6, *prob.globalBasis(), _p);
302
303
304
  prob.addVectorOperator(op6, _p);

  auto op7 = makeOperator(tag::partialtest_trial{0}, scalar);
305
  test(op7, *prob.globalBasis(), _p, _p);
306
307
308
  prob.addMatrixOperator(op7, _p, _p);

  auto op8 = makeOperator(tag::test_divtrialvec{}, scalar);
309
  test(op8, *prob.globalBasis(), _p, _v);
310
311
312
  prob.addMatrixOperator(op8, _p, _v);

  auto op9 = makeOperator(tag::test_gradtrial{}, vec);
313
  test(op9, *prob.globalBasis(), _p, _p);
314
315
316
  prob.addMatrixOperator(op9, _p, _p);

  auto op10 = makeOperator(tag::test_partialtrial{0}, scalar);
317
  test(op10, *prob.globalBasis(), _p, _p);
318
319
320
321
  prob.addMatrixOperator(op10, _p, _p);

  auto op11a = makeOperator(tag::testvec_gradtrial{}, scalar);
  // auto op11b = makeOperator(tag::testvec_gradtrial{}, mat); // TODO: needs to be implemented
322
  test(op11a, *prob.globalBasis(), _v, _p);
323
324
325
  prob.addMatrixOperator(op11a, _v, _p);

  auto op12 = makeOperator(tag::divtestvec_divtrialvec{}, scalar);
326
  test(op12, *prob.globalBasis(), _v, _v);
327
328
329
  prob.addMatrixOperator(op12, _v, _v);

  auto op13a = makeOperator(tag::gradtest_gradtrial{}, scalar);
330
  test(op13a, *prob.globalBasis(), _p, _p);
331
  prob.addMatrixOperator(op13a, _p, _p);
332
333
334

  auto op13b = makeOperator(tag::gradtest_gradtrial{}, mat);
  test(op13b, *prob.globalBasis(), _p, _p);
335
336
337
  prob.addMatrixOperator(op13b, _p, _p);

  auto op14 = makeOperator(tag::partialtest_partialtrial{0,0}, scalar);
338
  test(op14, *prob.globalBasis(), _p, _p);
339
340
341
  prob.addMatrixOperator(op14, _p, _p);

  auto op15 = makeOperator(tag::test{}, scalar);
342
  test(op15, *prob.globalBasis(), _p);
343
344
345
  prob.addVectorOperator(op15, _p);

  auto op16 = makeOperator(tag::test_trial{}, scalar);
346
  test(op16, *prob.globalBasis(), _p, _p);
347
348
349
  prob.addMatrixOperator(op16, _p, _p);

  auto op17 = makeOperator(tag::test_trialvec{}, vec);
350
  test(op17, *prob.globalBasis(), _p, _v);
351
352
353
  prob.addMatrixOperator(op17, _p, _v);

  auto op18 = makeOperator(tag::testvec{}, vec);
354
  test(op18, *prob.globalBasis(), _v);
355
356
357
  prob.addVectorOperator(op18, _v);

  auto op19 = makeOperator(tag::testvec_trial{}, vec);
358
  test(op19, *prob.globalBasis(), _v, _p);
359
360
361
  prob.addMatrixOperator(op19, _v, _p);

  auto op20a = makeOperator(tag::testvec_trialvec{}, scalar);
362
  test(op20a, *prob.globalBasis(), _v, _v);
363
  prob.addMatrixOperator(op20a, _v, _v);
364
365
366

  auto op20b = makeOperator(tag::testvec_trialvec{}, mat);
  test(op20b, *prob.globalBasis(), _v, _v);
367
368
369
370
371
  prob.addMatrixOperator(op20b, _v, _v);

  // some special operators

  auto opStokes = makeOperator(tag::stokes{}, scalar);
372
  test(opStokes, *prob.globalBasis(), makeTreePath(), makeTreePath());
373
374
375
376
  prob.addMatrixOperator(opStokes);
  prob.addMatrixOperator(opStokes, {}, {});

  auto opCDa = convectionDiffusion(scalar, scalar, scalar, scalar);
377
378
  test(opCDa, *prob.globalBasis(), _p, _p);
  test(opCDa, *prob.globalBasis(), _p);
379
380
381
382
  prob.addMatrixOperator(opCDa, _p, _p);
  prob.addVectorOperator(opCDa, _p);

  auto opCDb = convectionDiffusion(mat, vec, scalar, scalar, std::false_type{});
383
384
  test(opCDb, *prob.globalBasis(), _p, _p);
  test(opCDb, *prob.globalBasis(), _p);
385
386
387
  prob.addMatrixOperator(opCDb, _p, _p);
  prob.addVectorOperator(opCDb, _p);

388
389
390
391
392
  test_problem_operators(prob);

  auto prob2(prob);
  test_problem_operators(prob2);

393
394
  return 0;
}