#include <PreconditionerLowEnergy.h>
Inheritance diagram for Nektar::MultiRegions::PreconditionerLowEnergy:
Static Public Member Functions | |
| static PreconditionerSharedPtr | create (const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap) |
| Creates an instance of this class. More... | |
Static Public Attributes | |
| static std::string | className |
| Name of class. More... | |
Protected Member Functions | |
| virtual void | v_InitObject () override |
| virtual void | v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false) override |
| virtual void | v_BuildPreconditioner () override |
| Construct the low energy preconditioner from \(\mathbf{S}_{2}\). More... | |
| virtual DNekScalMatSharedPtr | v_TransformedSchurCompl (int n, int offset, const std::shared_ptr< DNekScalMat > &loc_mat) override |
| Set up the transformed block matrix system. More... | |
| virtual void | v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut) override |
| Transform the basis vector to low energy space. More... | |
| virtual void | v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut) override |
| transform the solution coeffiicents from low energy back to the original coefficient space. More... | |
| virtual void | v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override |
| Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis. More... | |
| virtual void | v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) override |
| Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space. More... | |
 Protected Member Functions inherited from Nektar::MultiRegions::Preconditioner | |
| virtual DNekScalMatSharedPtr | v_TransformedSchurCompl (int offset, int bndoffset, const std::shared_ptr< DNekScalMat > &loc_mat) |
| Get block elemental transposed transformation matrix \(\mathbf{R}^{T}\). More... | |
| virtual void | v_InitObject () |
| virtual void | v_DoPreconditioner (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const bool &isLocal=false) |
| Apply a preconditioner to the conjugate gradient method. More... | |
| virtual void | v_DoPreconditionerWithNonVertOutput (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput, const Array< OneD, NekDouble > &pNonVertOutput, Array< OneD, NekDouble > &pVertForce) |
| Apply a preconditioner to the conjugate gradient method with an output for non-vertex degrees of freedom. More... | |
| virtual void | v_DoTransformBasisToLowEnergy (Array< OneD, NekDouble > &pInOut) |
| Transform from original basis to low energy basis. More... | |
| virtual void | v_DoTransformCoeffsFromLowEnergy (Array< OneD, NekDouble > &pInOut) |
| Transform from low energy coeffs to orignal basis. More... | |
| virtual void | v_DoTransformCoeffsToLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) |
| Multiply by the block inverse transformation matrix. More... | |
| virtual void | v_DoTransformBasisFromLowEnergy (const Array< OneD, NekDouble > &pInput, Array< OneD, NekDouble > &pOutput) |
| Multiply by the block transposed inverse transformation matrix. More... | |
| virtual void | v_BuildPreconditioner () |
Protected Attributes | |
| DNekBlkMatSharedPtr | m_BlkMat |
| DNekBlkMatSharedPtr | m_RBlk |
| DNekBlkMatSharedPtr | m_InvRBlk |
| Array< OneD, NekDouble > | m_variablePmask |
| bool | m_signChange |
| std::vector< std::pair< int, int > > | m_sameBlock |
| Protected Attributes inherited from Nektar::MultiRegions::Preconditioner | |
| const std::weak_ptr< GlobalLinSys > | m_linsys |
| std::string | m_preconType |
| DNekMatSharedPtr | m_preconditioner |
| std::weak_ptr< AssemblyMap > | m_locToGloMap |
| LibUtilities::CommSharedPtr | m_comm |
Private Types | |
| typedef std::map< LibUtilities::ShapeType, DNekScalMatSharedPtr > | ShapeToDNekMap |
| typedef std::map< LibUtilities::ShapeType, LocalRegions::Expansion3DSharedPtr > | ShapeToExpMap |
| typedef std::map< LibUtilities::ShapeType, Array< OneD, unsigned int > > | ShapeToIntArrayMap |
| typedef std::map< LibUtilities::ShapeType, Array< OneD, Array< OneD, unsigned int > > > | ShapeToIntArrayArrayMap |
Detailed Description
This class implements low energy preconditioning for the conjugate gradient matrix solver.
Definition at line 53 of file PreconditionerLowEnergy.h.
Member Typedef Documentation
◆ ShapeToDNekMap
|
private |
Definition at line 123 of file PreconditionerLowEnergy.h.
◆ ShapeToExpMap
|
private |
Definition at line 126 of file PreconditionerLowEnergy.h.
◆ ShapeToIntArrayArrayMap
|
private |
Definition at line 131 of file PreconditionerLowEnergy.h.
◆ ShapeToIntArrayMap
|
private |
Definition at line 128 of file PreconditionerLowEnergy.h.
Constructor & Destructor Documentation
◆ PreconditionerLowEnergy()
| Nektar::MultiRegions::PreconditionerLowEnergy::PreconditionerLowEnergy | ( | const std::shared_ptr< GlobalLinSys > & | plinsys, |
| const AssemblyMapSharedPtr & | pLocToGloMap | ||
| ) |
Definition at line 65 of file PreconditionerLowEnergy.cpp.
68 : Preconditioner(plinsys, pLocToGloMap)
69{
70}
Preconditioner(const std::shared_ptr< GlobalLinSys > &plinsys, const AssemblyMapSharedPtr &pLocToGloMap)
Definition: Preconditioner.cpp:58
◆ ~PreconditionerLowEnergy()
|
inlinevirtual |
Definition at line 76 of file PreconditionerLowEnergy.h.
77 {
78 }
Member Function Documentation
◆ create()
|
inlinestatic |
Creates an instance of this class.
Definition at line 57 of file PreconditionerLowEnergy.h.
60 {
63 plinsys, pLocToGloMap);
64 p->InitObject();
66 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168
std::shared_ptr< Preconditioner > PreconditionerSharedPtr
Definition: GlobalLinSys.h:60
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.
◆ CreateRefHexGeom()
|
private |
Sets up the reference hexahedral element needed to construct a low energy basis.
Definition at line 1542 of file PreconditionerLowEnergy.cpp.
1543{
1544 ////////////////////////////////
1545 // Set up Hexahedron vertices //
1546 ////////////////////////////////
1547
1548 const int three = 3;
1549
1550 const int nVerts = 8;
1551 const double point[][3] = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0},
1552 {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}};
1553
1554 // Populate the list of verts
1555 SpatialDomains::PointGeomSharedPtr verts[8];
1556 for (int i = 0; i < nVerts; ++i)
1557 {
1558 verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
1559 three, i, point[i][0], point[i][1], point[i][2]);
1560 }
1561
1562 /////////////////////////////
1563 // Set up Hexahedron Edges //
1564 /////////////////////////////
1565
1566 // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1567 const int nEdges = 12;
1568 const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {2, 3}, {0, 3},
1569 {0, 4}, {1, 5}, {2, 6}, {3, 7},
1570 {4, 5}, {5, 6}, {6, 7}, {4, 7}};
1571
1572 // Populate the list of edges
1573 SpatialDomains::SegGeomSharedPtr edges[nEdges];
1574 for (int i = 0; i < nEdges; ++i)
1575 {
1576 SpatialDomains::PointGeomSharedPtr vertsArray[2];
1577 for (int j = 0; j < 2; ++j)
1578 {
1579 vertsArray[j] = verts[vertexConnectivity[i][j]];
1580 }
1581 edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
1582 i, three, vertsArray);
1583 }
1584
1585 /////////////////////////////
1586 // Set up Hexahedron faces //
1587 /////////////////////////////
1588
1589 const int nFaces = 6;
1590 const int edgeConnectivity[][4] = {{0, 1, 2, 3}, {0, 5, 8, 4},
1591 {1, 6, 9, 5}, {2, 7, 10, 6},
1592 {3, 7, 11, 4}, {8, 9, 10, 11}};
1593
1594 // Populate the list of faces
1595 SpatialDomains::QuadGeomSharedPtr faces[nFaces];
1596 for (int i = 0; i < nFaces; ++i)
1597 {
1598 SpatialDomains::SegGeomSharedPtr edgeArray[4];
1599 for (int j = 0; j < 4; ++j)
1600 {
1601 edgeArray[j] = edges[edgeConnectivity[i][j]];
1602 }
1603 faces[i] = MemoryManager<SpatialDomains::QuadGeom>::AllocateSharedPtr(
1604 i, edgeArray);
1605 }
1606
1607 SpatialDomains::HexGeomSharedPtr geom =
1608 MemoryManager<SpatialDomains::HexGeom>::AllocateSharedPtr(0, faces);
1609
1610 return geom;
1611}
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:47
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:87
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
References Nektar::MemoryManager< DataType >::AllocateSharedPtr().
Referenced by SetUpReferenceElements().
◆ CreateRefPrismGeom()
|
private |
Sets up the reference prismatic element needed to construct a low energy basis.
Definition at line 1279 of file PreconditionerLowEnergy.cpp.
1280{
1281 //////////////////////////
1282 // Set up Prism element //
1283 //////////////////////////
1284
1285 const int three = 3;
1286 const int nVerts = 6;
1287 const double point[][3] = {
1288 {-1, -1, 0},
1289 {1, -1, 0},
1290 {1, 1, 0},
1291 {-1, 1, 0},
1294 };
1295
1296 // std::shared_ptr<SpatialDomains::PointGeom> verts[6];
1297 SpatialDomains::PointGeomSharedPtr verts[6];
1298 for (int i = 0; i < nVerts; ++i)
1299 {
1300 verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
1301 three, i, point[i][0], point[i][1], point[i][2]);
1302 }
1303 const int nEdges = 9;
1304 const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {3, 2}, {0, 3}, {0, 4},
1305 {1, 4}, {2, 5}, {3, 5}, {4, 5}};
1306
1307 // Populate the list of edges
1308 SpatialDomains::SegGeomSharedPtr edges[nEdges];
1309 for (int i = 0; i < nEdges; ++i)
1310 {
1311 SpatialDomains::PointGeomSharedPtr vertsArray[2];
1312 for (int j = 0; j < 2; ++j)
1313 {
1314 vertsArray[j] = verts[vertexConnectivity[i][j]];
1315 }
1316 edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
1317 i, three, vertsArray);
1318 }
1319
1320 ////////////////////////
1321 // Set up Prism faces //
1322 ////////////////////////
1323
1324 const int nFaces = 5;
1325 // quad-edge connectivity base-face0, vertical-quadface2, vertical-quadface4
1326 const int quadEdgeConnectivity[][4] = {
1327 {0, 1, 2, 3}, {1, 6, 8, 5}, {3, 7, 8, 4}};
1328 // QuadId ordered as 0, 1, 2, otherwise return false
1329 const int quadId[] = {0, -1, 1, -1, 2};
1330
1331 // triangle-edge connectivity side-triface-1, side triface-3
1332 const int triEdgeConnectivity[][3] = {{0, 5, 4}, {2, 6, 7}};
1333 // TriId ordered as 0, 1, otherwise return false
1334 const int triId[] = {-1, 0, -1, 1, -1};
1335
1336 // Populate the list of faces
1337 SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1338 for (int f = 0; f < nFaces; ++f)
1339 {
1340 if (f == 1 || f == 3)
1341 {
1342 int i = triId[f];
1343 SpatialDomains::SegGeomSharedPtr edgeArray[3];
1344 for (int j = 0; j < 3; ++j)
1345 {
1346 edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1347 }
1348 faces[f] =
1350 f, edgeArray);
1351 }
1352 else
1353 {
1354 int i = quadId[f];
1355 SpatialDomains::SegGeomSharedPtr edgeArray[4];
1356 for (int j = 0; j < 4; ++j)
1357 {
1358 edgeArray[j] = edges[quadEdgeConnectivity[i][j]];
1359 }
1360 faces[f] =
1362 f, edgeArray);
1363 }
1364 }
1365
1366 SpatialDomains::PrismGeomSharedPtr geom =
1367 MemoryManager<SpatialDomains::PrismGeom>::AllocateSharedPtr(0, faces);
1368
1369 return geom;
1370}
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:85
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and tinysimd::sqrt().
Referenced by SetUpReferenceElements().
◆ CreateRefPyrGeom()
|
private |
Sets up the reference prismatic element needed to construct a low energy basis mapping arrays.
Definition at line 1376 of file PreconditionerLowEnergy.cpp.
1377{
1378 //////////////////////////
1379 // Set up Pyramid element //
1380 //////////////////////////
1381
1382 const int nVerts = 5;
1383 const double point[][3] = {{-1, -1, 0},
1384 {1, -1, 0},
1385 {1, 1, 0},
1386 {-1, 1, 0},
1388
1389 // boost::shared_ptr<SpatialDomains::PointGeom> verts[6];
1390 const int three = 3;
1391 SpatialDomains::PointGeomSharedPtr verts[5];
1392 for (int i = 0; i < nVerts; ++i)
1393 {
1394 verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
1395 three, i, point[i][0], point[i][1], point[i][2]);
1396 }
1397 const int nEdges = 8;
1398 const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0},
1399 {0, 4}, {1, 4}, {2, 4}, {3, 4}};
1400
1401 // Populate the list of edges
1402 SpatialDomains::SegGeomSharedPtr edges[nEdges];
1403 for (int i = 0; i < nEdges; ++i)
1404 {
1405 SpatialDomains::PointGeomSharedPtr vertsArray[2];
1406 for (int j = 0; j < 2; ++j)
1407 {
1408 vertsArray[j] = verts[vertexConnectivity[i][j]];
1409 }
1410 edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
1411 i, three, vertsArray);
1412 }
1413
1414 ////////////////////////
1415 // Set up Pyramid faces //
1416 ////////////////////////
1417
1418 const int nFaces = 5;
1419 // quad-edge connectivity base-face0,
1420 const int quadEdgeConnectivity[][4] = {{0, 1, 2, 3}};
1421
1422 // triangle-edge connectivity side-triface-1, side triface-2
1423 const int triEdgeConnectivity[][3] = {
1424 {0, 5, 4}, {1, 6, 5}, {2, 7, 6}, {3, 4, 7}};
1425
1426 // Populate the list of faces
1427 SpatialDomains::Geometry2DSharedPtr faces[nFaces];
1428 for (int f = 0; f < nFaces; ++f)
1429 {
1430 if (f == 0)
1431 {
1432 SpatialDomains::SegGeomSharedPtr edgeArray[4];
1433 for (int j = 0; j < 4; ++j)
1434 {
1435 edgeArray[j] = edges[quadEdgeConnectivity[f][j]];
1436 }
1437
1438 faces[f] =
1440 f, edgeArray);
1441 }
1442 else
1443 {
1444 int i = f - 1;
1445 SpatialDomains::SegGeomSharedPtr edgeArray[3];
1446 for (int j = 0; j < 3; ++j)
1447 {
1448 edgeArray[j] = edges[triEdgeConnectivity[i][j]];
1449 }
1450 faces[f] =
1452 f, edgeArray);
1453 }
1454 }
1455
1456 SpatialDomains::PyrGeomSharedPtr geom =
1457 MemoryManager<SpatialDomains::PyrGeom>::AllocateSharedPtr(0, faces);
1458
1459 return geom;
1460}
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:77
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and tinysimd::sqrt().
Referenced by SetUpReferenceElements().
◆ CreateRefTetGeom()
|
private |
Sets up the reference tretrahedral element needed to construct a low energy basis.
Definition at line 1466 of file PreconditionerLowEnergy.cpp.
1467{
1468 /////////////////////////////////
1469 // Set up Tetrahedron vertices //
1470 /////////////////////////////////
1471
1472 int i, j;
1473 const int three = 3;
1474 const int nVerts = 4;
1479
1480 std::shared_ptr<SpatialDomains::PointGeom> verts[4];
1481 for (i = 0; i < nVerts; ++i)
1482 {
1483 verts[i] = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
1484 three, i, point[i][0], point[i][1], point[i][2]);
1485 }
1486
1487 //////////////////////////////
1488 // Set up Tetrahedron Edges //
1489 //////////////////////////////
1490
1491 // SegGeom (int id, const int coordim), EdgeComponent(id, coordim)
1492 const int nEdges = 6;
1493 const int vertexConnectivity[][2] = {{0, 1}, {1, 2}, {0, 2},
1494 {0, 3}, {1, 3}, {2, 3}};
1495
1496 // Populate the list of edges
1497 SpatialDomains::SegGeomSharedPtr edges[nEdges];
1498 for (i = 0; i < nEdges; ++i)
1499 {
1500 std::shared_ptr<SpatialDomains::PointGeom> vertsArray[2];
1501 for (j = 0; j < 2; ++j)
1502 {
1503 vertsArray[j] = verts[vertexConnectivity[i][j]];
1504 }
1505
1506 edges[i] = MemoryManager<SpatialDomains::SegGeom>::AllocateSharedPtr(
1507 i, three, vertsArray);
1508 }
1509
1510 //////////////////////////////
1511 // Set up Tetrahedron faces //
1512 //////////////////////////////
1513
1514 const int nFaces = 4;
1515 const int edgeConnectivity[][3] = {
1516 {0, 1, 2}, {0, 4, 3}, {1, 5, 4}, {2, 5, 3}};
1517
1518 // Populate the list of faces
1519 SpatialDomains::TriGeomSharedPtr faces[nFaces];
1520 for (i = 0; i < nFaces; ++i)
1521 {
1522 SpatialDomains::SegGeomSharedPtr edgeArray[3];
1523 for (j = 0; j < 3; ++j)
1524 {
1525 edgeArray[j] = edges[edgeConnectivity[i][j]];
1526 }
1527
1528 faces[i] = MemoryManager<SpatialDomains::TriGeom>::AllocateSharedPtr(
1529 i, edgeArray);
1530 }
1531
1532 SpatialDomains::TetGeomSharedPtr geom =
1533 MemoryManager<SpatialDomains::TetGeom>::AllocateSharedPtr(0, faces);
1534
1535 return geom;
1536}
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:85
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and tinysimd::sqrt().
Referenced by SetUpReferenceElements().
◆ CreateVariablePMask()
Create the mask array
Definition at line 1252 of file PreconditionerLowEnergy.cpp.
1253{
1255 unsigned int nLocBnd = asmMap->GetNumLocalBndCoeffs();
1256
1257 m_signChange = asmMap->GetSignChange();
1258
1259 // Construct a mask array
1260 m_variablePmask = Array<OneD, NekDouble>(nLocBnd, 1.0);
1261
1263 {
1264 unsigned int i;
1266 asmMap->GetLocalToGlobalBndSign();
1267
1268 for (i = 0; i < nLocBnd; ++i)
1269 {
1271 }
1272 }
1273}
std::weak_ptr< AssemblyMap > m_locToGloMap
Definition: Preconditioner.h:125
Array< OneD, NekDouble > m_variablePmask
Definition: PreconditionerLowEnergy.h:85
bool m_signChange
Definition: PreconditionerLowEnergy.h:87
References Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_signChange, m_variablePmask, and sign.
Referenced by v_InitObject().
◆ ExtractLocMat()
|
private |
Definition at line 2265 of file PreconditionerLowEnergy.cpp.
2269{
2270 NekDouble val;
2271 NekDouble zero = 0.0;
2272
2273 int nRows = locExp->NumBndryCoeffs();
2274 DNekMatSharedPtr newmat =
2276
2277 Array<OneD, unsigned int> vlocmap;
2278 Array<OneD, Array<OneD, unsigned int>> elocmap;
2279 Array<OneD, Array<OneD, unsigned int>> flocmap;
2280
2281 locExp->GetInverseBoundaryMaps(vlocmap, elocmap, flocmap);
2282
2283 // fill diagonal
2284 for (int i = 0; i < nRows; ++i)
2285 {
2286 val = 1.0;
2287 newmat->SetValue(i, i, val);
2288 }
2289
2290 int nverts = locExp->GetNverts();
2291 int nedges = locExp->GetNedges();
2292 int nfaces = locExp->GetNtraces();
2293
2294 // fill vertex to edge coupling
2295 for (int e = 0; e < nedges; ++e)
2296 {
2297 int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
2298
2299 for (int v = 0; v < nverts; ++v)
2300 {
2301 for (int i = 0; i < nEdgeInteriorCoeffs; ++i)
2302 {
2303 val = (*maxRmat)(vmap[v], emap[e][i]);
2304 newmat->SetValue(vlocmap[v], elocmap[e][i], val);
2305 }
2306 }
2307 }
2308
2309 for (int f = 0; f < nfaces; ++f)
2310 {
2311 // Get details to extrac this face from max reference matrix
2312 StdRegions::Orientation FwdOrient =
2314 int m0, m1; // Local face expansion orders.
2315
2316 int nFaceInteriorCoeffs = locExp->GetTraceIntNcoeffs(f);
2317
2318 locExp->GetTraceNumModes(f, m0, m1, FwdOrient);
2319
2320 Array<OneD, unsigned int> fmapRmat =
2321 maxExp->GetTraceInverseBoundaryMap(f, FwdOrient, m0, m1);
2322
2323 // fill in vertex to face coupling
2324 for (int v = 0; v < nverts; ++v)
2325 {
2326 for (int i = 0; i < nFaceInteriorCoeffs; ++i)
2327 {
2328 val = (*maxRmat)(vmap[v], fmapRmat[i]);
2329 newmat->SetValue(vlocmap[v], flocmap[f][i], val);
2330 }
2331 }
2332
2333 // fill in edges to face coupling
2334 for (int e = 0; e < nedges; ++e)
2335 {
2336 int nEdgeInteriorCoeffs = locExp->GetEdgeNcoeffs(e) - 2;
2337
2338 for (int j = 0; j < nEdgeInteriorCoeffs; ++j)
2339 {
2340 for (int i = 0; i < nFaceInteriorCoeffs; ++i)
2341 {
2342 val = (*maxRmat)(emap[e][j], fmapRmat[i]);
2343 newmat->SetValue(elocmap[e][j], flocmap[f][i], val);
2344 }
2345 }
2346 }
2347 }
2348
2349 return newmat;
2350}
@ eDir1FwdDir1_Dir2FwdDir2
Definition: StdRegions.hpp:422
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::eDir1FwdDir1_Dir2FwdDir2, and Nektar::eFULL.
Referenced by SetupBlockTransformationMatrix().
◆ ReSetPrismMaxRMat()
|
private |
Definition at line 2049 of file PreconditionerLowEnergy.cpp.
2056{
2057 int nRows = PrismExp->NumBndryCoeffs();
2058 NekDouble val;
2059 NekDouble zero = 0.0;
2060 DNekMatSharedPtr newmat =
2062
2063 int nfacemodes;
2064
2065 if (UseTetOnly)
2066 {
2067 // copy existing system
2068 for (int i = 0; i < nRows; ++i)
2069 {
2070 for (int j = 0; j < nRows; ++j)
2071 {
2072 val = (*maxRmat[ePrism])(i, j);
2073 newmat->SetValue(i, j, val);
2074 }
2075 }
2076
2077 // Reset vertex to edge mapping from tet.
2078 const int VETetVert[][2] = {{0, 0}, {1, 1}, {1, 1},
2079 {0, 0}, {3, 3}, {3, 3}};
2080 const int VETetEdge[][2] = {{0, 3}, {0, 4}, {0, 4},
2081 {0, 3}, {3, 4}, {4, 3}};
2082 const int VEPrismEdge[][2] = {{0, 4}, {0, 5}, {2, 6},
2083 {2, 7}, {4, 5}, {6, 7}};
2084
2085 // fill vertex to edge coupling
2086 for (int v = 0; v < 6; ++v)
2087 {
2088 for (int e = 0; e < 2; ++e)
2089 {
2090 for (int i = 0; i < nummodesmax - 2; ++i)
2091 {
2092 // Note this is using wrong shape but gives
2093 // answer that seems to have correct error!
2094 val = (*maxRmat[eTetrahedron])(
2095 vertMapMaxR[eTetrahedron][VETetVert[v][e]],
2096 edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i]);
2097 newmat->SetValue(vertMapMaxR[ePrism][v],
2098 edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
2099 val);
2100 }
2101 }
2102 }
2103 }
2104 else
2105 {
2106 // set diagonal to 1
2107 for (int i = 0; i < nRows; ++i)
2108 {
2109 newmat->SetValue(i, i, 1.0);
2110 }
2111
2112 // Set vertex to edge mapping from Hex.
2113
2114 // The following lists specify the number of adjacent
2115 // edges to each vertex (nadj) then the Hex vert to
2116 // use for each prism ver in the vert-edge map (VEVert)
2117 // followed by the hex edge to use for each prism edge
2118 // in the vert-edge map (VEEdge)
2119 const int VEHexVert[][3] = {{0, 0, 0}, {1, 1, 1}, {2, 2, 2},
2120 {3, 3, 3}, {4, 5, 5}, {6, 7, 7}};
2121 const int VEHexEdge[][3] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6},
2122 {2, 3, 7}, {4, 5, 9}, {6, 7, 11}};
2123 const int VEPrismEdge[][3] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6},
2124 {2, 3, 7}, {4, 5, 8}, {6, 7, 8}};
2125
2126 // fill vertex to edge coupling
2127 for (int v = 0; v < 6; ++v)
2128 {
2129 for (int e = 0; e < 3; ++e)
2130 {
2131 for (int i = 0; i < nummodesmax - 2; ++i)
2132 {
2133 // Note this is using wrong shape but gives
2134 // answer that seems to have correct error!
2135 val = (*maxRmat[eHexahedron])(
2136 vertMapMaxR[eHexahedron][VEHexVert[v][e]],
2137 edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
2138 newmat->SetValue(vertMapMaxR[ePrism][v],
2139 edgeMapMaxR[ePrism][VEPrismEdge[v][e]][i],
2140 val);
2141 }
2142 }
2143 }
2144
2145 // Setup vertex to face mapping from Hex
2146 const int VFHexVert[][2] = {{0, 0}, {1, 1}, {4, 5},
2147 {2, 2}, {3, 3}, {6, 7}};
2148 const int VFHexFace[][2] = {{0, 4}, {0, 2}, {4, 2},
2149 {0, 2}, {0, 4}, {2, 4}};
2150
2151 const int VQFPrismVert[][2] = {{0, 0}, {1, 1}, {4, 4},
2152 {2, 2}, {3, 3}, {5, 5}};
2153 const int VQFPrismFace[][2] = {{0, 4}, {0, 2}, {4, 2},
2154 {0, 2}, {0, 4}, {2, 4}};
2155
2156 nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
2157 // Replace two Quad faces on every vertex
2158 for (int v = 0; v < 6; ++v)
2159 {
2160 for (int f = 0; f < 2; ++f)
2161 {
2162 for (int i = 0; i < nfacemodes; ++i)
2163 {
2164 val = (*maxRmat[eHexahedron])(
2165 vertMapMaxR[eHexahedron][VFHexVert[v][f]],
2166 faceMapMaxR[eHexahedron][VFHexFace[v][f]][i]);
2167 newmat->SetValue(vertMapMaxR[ePrism][VQFPrismVert[v][f]],
2168 faceMapMaxR[ePrism][VQFPrismFace[v][f]][i],
2169 val);
2170 }
2171 }
2172 }
2173
2174 // Mapping of Hex Edge-Face mappings to Prism Edge-Face Mappings
2175 const int nadjface[] = {1, 2, 1, 2, 1, 1, 1, 1, 2};
2176 const int EFHexEdge[][2] = {{0, -1}, {1, 1}, {2, -1}, {3, 3}, {4, -1},
2177 {5, -1}, {6, -1}, {7, -1}, {9, 11}};
2178 const int EFHexFace[][2] = {{0, -1}, {0, 2}, {0, -1}, {0, 4}, {4, -1},
2179 {2, -1}, {2, -1}, {4, -1}, {2, 4}};
2180 const int EQFPrismEdge[][2] = {{0, -1}, {1, 1}, {2, -1},
2181 {3, 3}, {4, -1}, {5, -1},
2182 {6, -1}, {7, -1}, {8, 8}};
2183 const int EQFPrismFace[][2] = {{0, -1}, {0, 2}, {0, -1},
2184 {0, 4}, {4, -1}, {2, -1},
2185 {2, -1}, {4, -1}, {2, 4}};
2186
2187 // all base edges are coupled to face 0
2188 nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
2189 for (int e = 0; e < 9; ++e)
2190 {
2191 for (int f = 0; f < nadjface[e]; ++f)
2192 {
2193 for (int i = 0; i < nummodesmax - 2; ++i)
2194 {
2195 for (int j = 0; j < nfacemodes; ++j)
2196 {
2197 int edgemapid =
2198 edgeMapMaxR[eHexahedron][EFHexEdge[e][f]][i];
2199 int facemapid =
2200 faceMapMaxR[eHexahedron][EFHexFace[e][f]][j];
2201
2202 val = (*maxRmat[eHexahedron])(edgemapid, facemapid);
2203
2204 int edgemapid1 =
2205 edgeMapMaxR[ePrism][EQFPrismEdge[e][f]][i];
2206 int facemapid1 =
2207 faceMapMaxR[ePrism][EQFPrismFace[e][f]][j];
2208 newmat->SetValue(edgemapid1, facemapid1, val);
2209 }
2210 }
2211 }
2212 }
2213 }
2214
2215 const int VFTetVert[] = {0, 1, 3, 1, 0, 3};
2216 const int VFTetFace[] = {1, 1, 1, 1, 1, 1};
2217 const int VTFPrismVert[] = {0, 1, 4, 2, 3, 5};
2218 const int VTFPrismFace[] = {1, 1, 1, 3, 3, 3};
2219
2220 // Handle all triangular faces from tetrahedron
2221 nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
2222 for (int v = 0; v < 6; ++v)
2223 {
2224 for (int i = 0; i < nfacemodes; ++i)
2225 {
2226 val = (*maxRmat[eTetrahedron])(
2227 vertMapMaxR[eTetrahedron][VFTetVert[v]],
2228 faceMapMaxR[eTetrahedron][VFTetFace[v]][i]);
2229
2230 newmat->SetValue(vertMapMaxR[ePrism][VTFPrismVert[v]],
2231 faceMapMaxR[ePrism][VTFPrismFace[v]][i], val);
2232 }
2233 }
2234
2235 // Mapping of Tet Edge-Face mappings to Prism Edge-Face Mappings
2236 const int EFTetEdge[] = {0, 3, 4, 0, 4, 3};
2237 const int EFTetFace[] = {1, 1, 1, 1, 1, 1};
2238 const int ETFPrismEdge[] = {0, 4, 5, 2, 6, 7};
2239 const int ETFPrismFace[] = {1, 1, 1, 3, 3, 3};
2240
2241 // handle all edge to triangular faces from tetrahedron
2242 // (only 6 this time)
2243 nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
2244 for (int e = 0; e < 6; ++e)
2245 {
2246 for (int i = 0; i < nummodesmax - 2; ++i)
2247 {
2248 for (int j = 0; j < nfacemodes; ++j)
2249 {
2252 val = (*maxRmat[eTetrahedron])(edgemapid, facemapid);
2253
2254 newmat->SetValue(edgeMapMaxR[ePrism][ETFPrismEdge[e]][i],
2255 faceMapMaxR[ePrism][ETFPrismFace[e]][j], val);
2256 }
2257 }
2258 }
2259
2260 DNekScalMatSharedPtr PrismR =
2261 MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
2262 maxRmat[ePrism] = PrismR;
2263}
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:68
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePrism, and Nektar::LibUtilities::eTetrahedron.
Referenced by SetUpReferenceElements().
◆ ReSetTetMaxRMat()
|
private |
Definition at line 1990 of file PreconditionerLowEnergy.cpp.
1996{
1997 boost::ignore_unused(faceMapMaxR);
1998
1999 int nRows = TetExp->NumBndryCoeffs();
2000 NekDouble val;
2001 NekDouble zero = 0.0;
2002 DNekMatSharedPtr newmat =
2004
2005 // copy existing system
2006 for (int i = 0; i < nRows; ++i)
2007 {
2008 for (int j = 0; j < nRows; ++j)
2009 {
2010 val = (*maxRmat[eTetrahedron])(i, j);
2011 newmat->SetValue(i, j, val);
2012 }
2013 }
2014
2015 // The following lists specify the number of adjacent
2016 // edges to each vertex (nadj) then the Hex vert to
2017 // use for each pyramid ver in the vert-edge map (VEVert)
2018 // followed by the hex edge to use for each Tet edge
2019 // in the vert-edge map (VEEdge)
2020 const int VEHexVert[][4] = {{0, 0, 0}, {1, 1, 1}, {2, 2, 2}, {4, 5, 6}};
2021 const int VEHexEdge[][4] = {{0, 3, 4}, {0, 1, 5}, {1, 2, 6}, {4, 5, 6}};
2022 const int VETetEdge[][4] = {{0, 2, 3}, {0, 1, 4}, {1, 2, 5}, {3, 4, 5}};
2023
2024 // fill vertex to edge coupling
2025 for (int v = 0; v < 4; ++v)
2026 {
2027 for (int e = 0; e < 3; ++e)
2028 {
2029 for (int i = 0; i < nummodesmax - 2; ++i)
2030 {
2031 // Note this is using wrong shape but gives
2032 // answer that seems to have correct error!
2033 val = (*maxRmat[eHexahedron])(
2034 vertMapMaxR[eHexahedron][VEHexVert[v][e]],
2035 edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
2036 newmat->SetValue(vertMapMaxR[eTetrahedron][v],
2037 edgeMapMaxR[eTetrahedron][VETetEdge[v][e]][i],
2038 val);
2039 }
2040 }
2041 }
2042
2043 DNekScalMatSharedPtr TetR =
2044 MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
2045
2046 maxRmat[eTetrahedron] = TetR;
2047}
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, and Nektar::LibUtilities::eTetrahedron.
Referenced by SetUpReferenceElements().
◆ SetupBlockTransformationMatrix()
Set a block transformation matrices for each element type. These are needed in routines that transform the schur complement matrix to and from the low energy basis.
Definition at line 936 of file PreconditionerLowEnergy.cpp.
937{
938 std::shared_ptr<MultiRegions::ExpList> expList =
939 ((m_linsys.lock())->GetLocMat()).lock();
940 LocalRegions::Expansion3DSharedPtr locExp;
941 LocalRegions::Expansion3DSharedPtr locExpSav;
942 map<int, int> EdgeSize;
943
944 int n;
945
946 std::map<ShapeType, DNekScalMatSharedPtr> maxRmat;
947 map<ShapeType, LocalRegions::Expansion3DSharedPtr> maxElmt;
948 map<ShapeType, Array<OneD, unsigned int>> vertMapMaxR;
949 map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> edgeMapMaxR;
950
951 // Sets up reference element and builds transformation matrix for
952 // maximum polynomial order meshes
953 SetUpReferenceElements(maxRmat, maxElmt, vertMapMaxR, edgeMapMaxR);
954
955 const Array<OneD, const unsigned int> &nbdry_size =
956 m_locToGloMap.lock()->GetNumLocalBndCoeffsPerPatch();
957
958 int n_exp = expList->GetNumElmts();
959
961
962 // Variants of R matrices required for low energy preconditioning
964 nbdry_size, nbdry_size, blkmatStorage);
966 nbdry_size, nbdry_size, blkmatStorage);
967
968 DNekMatSharedPtr rmat, invrmat;
969
970 int offset = 0;
971
972 // Set up transformation matrices whilst checking to see if
973 // consecutive matrices are the same and if so reuse the
974 // matrices and store how many consecutive offsets there
975 // are
976 for (n = 0; n < n_exp; ++n)
977 {
978 locExp = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
979
980 ShapeType eltype = locExp->DetShapeType();
981
982 int nbndcoeffs = locExp->NumBndryCoeffs();
983
985 {
986 rmat = ExtractLocMat(locExp, maxRmat[eltype], maxElmt[eltype],
987 vertMapMaxR[eltype], edgeMapMaxR[eltype]);
988 // Block R matrix
989 m_RBlk->SetBlock(n, n, rmat);
990
991 invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
992 invrmat->Invert();
993
994 // Block inverse R matrix
995 m_InvRBlk->SetBlock(n, n, invrmat);
996
997 m_sameBlock.push_back(pair<int, int>(1, nbndcoeffs));
998 locExpSav = locExp;
999 }
1000 else
1001 {
1002 bool reuse = true;
1003
1004 // check to see if same as previous matrix and
1005 // reuse if we can
1006 for (int i = 0; i < 3; ++i)
1007 {
1008 if (locExpSav->GetBasis(i) != locExp->GetBasis(i))
1009 {
1010 reuse = false;
1011 break;
1012 }
1013 }
1014
1015 if (reuse)
1016 {
1017 m_RBlk->SetBlock(n, n, rmat);
1018 m_InvRBlk->SetBlock(n, n, invrmat);
1019
1020 m_sameBlock[offset] =
1021 (pair<int, int>(m_sameBlock[offset].first + 1, nbndcoeffs));
1022 }
1023 else
1024 {
1025 rmat = ExtractLocMat(locExp, maxRmat[eltype], maxElmt[eltype],
1026 vertMapMaxR[eltype], edgeMapMaxR[eltype]);
1027
1028 // Block R matrix
1029 m_RBlk->SetBlock(n, n, rmat);
1030
1031 invrmat = MemoryManager<DNekMat>::AllocateSharedPtr(*rmat);
1032 invrmat->Invert();
1033 // Block inverse R matrix
1034 m_InvRBlk->SetBlock(n, n, invrmat);
1035
1036 m_sameBlock.push_back(pair<int, int>(1, nbndcoeffs));
1037 offset++;
1038 locExpSav = locExp;
1039 }
1040 }
1041 }
1042}
const std::weak_ptr< GlobalLinSys > m_linsys
Definition: Preconditioner.h:122
std::vector< std::pair< int, int > > m_sameBlock
Definition: PreconditionerLowEnergy.h:91
DNekBlkMatSharedPtr m_InvRBlk
Definition: PreconditionerLowEnergy.h:83
DNekBlkMatSharedPtr m_RBlk
Definition: PreconditionerLowEnergy.h:82
DNekMatSharedPtr ExtractLocMat(LocalRegions::Expansion3DSharedPtr &locExp, DNekScalMatSharedPtr &maxRmat, LocalRegions::Expansion3DSharedPtr &expMax, Array< OneD, unsigned int > &vertMapMaxR, Array< OneD, Array< OneD, unsigned int > > &edgeMapMaxR)
Definition: PreconditionerLowEnergy.cpp:2265
void SetUpReferenceElements(ShapeToDNekMap &maxRmat, ShapeToExpMap &maxElmt, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR)
Loop expansion and determine different variants of the transformation matrix.
Definition: PreconditionerLowEnergy.cpp:1619
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:51
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::StdRegions::StdExpansion::DetShapeType(), Nektar::eDIAGONAL, ExtractLocMat(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, and SetUpReferenceElements().
Referenced by v_InitObject().
◆ SetUpPyrMaxRMat()
|
private |
Definition at line 1825 of file PreconditionerLowEnergy.cpp.
1831{
1832 int nRows = PyrExp->NumBndryCoeffs();
1833 NekDouble val;
1834 NekDouble zero = 0.0;
1835 DNekMatSharedPtr newmat =
1837
1838 // set diagonal to 1
1839 for (int i = 0; i < nRows; ++i)
1840 {
1841 newmat->SetValue(i, i, 1.0);
1842 }
1843
1844 // The following lists specify the number of adjacent
1845 // edges to each vertex (nadj) then the Hex vert to
1846 // use for each pyramid ver in the vert-edge map (VEVert)
1847 // followed by the hex edge to use for each pyramid edge
1848 // in the vert-edge map (VEEdge)
1849 const int nadjedge[] = {3, 3, 3, 3, 4};
1850 const int VEHexVert[][4] = {{0, 0, 0, -1},
1851 {1, 1, 1, -1},
1852 {2, 2, 2, -1},
1853 {3, 3, 3, -1},
1854 {4, 5, 6, 7}};
1855 const int VEHexEdge[][4] = {{0, 3, 4, -1},
1856 {0, 1, 5, -1},
1857 {1, 2, 6, -1},
1858 {2, 3, 7, -1},
1859 {4, 5, 6, 7}};
1860 const int VEPyrEdge[][4] = {{0, 3, 4, -1},
1861 {0, 1, 5, -1},
1862 {1, 2, 6, -1},
1863 {2, 3, 7, -1},
1864 {4, 5, 6, 7}};
1865
1866 // fill vertex to edge coupling
1867 for (int v = 0; v < 5; ++v)
1868 {
1869 for (int e = 0; e < nadjedge[v]; ++e)
1870 {
1871 for (int i = 0; i < nummodesmax - 2; ++i)
1872 {
1873 // Note this is using wrong shape but gives
1874 // answer that seems to have correct error!
1875 val = (*maxRmat[eHexahedron])(
1876 vertMapMaxR[eHexahedron][VEHexVert[v][e]],
1877 edgeMapMaxR[eHexahedron][VEHexEdge[v][e]][i]);
1878 newmat->SetValue(vertMapMaxR[ePyramid][v],
1879 edgeMapMaxR[ePyramid][VEPyrEdge[v][e]][i],
1880 val);
1881 }
1882 }
1883 }
1884
1885 int nfacemodes;
1886 nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
1887 // First four verties are all adjacent to base face
1888 for (int v = 0; v < 4; ++v)
1889 {
1890 for (int i = 0; i < nfacemodes; ++i)
1891 {
1893 faceMapMaxR[eHexahedron][0][i]);
1894 newmat->SetValue(vertMapMaxR[ePyramid][v],
1895 faceMapMaxR[ePyramid][0][i], val);
1896 }
1897 }
1898
1899 const int nadjface[] = {2, 2, 2, 2, 4};
1900 const int VFTetVert[][4] = {{0, 0, -1, -1},
1901 {1, 1, -1, -1},
1902 {2, 2, -1, -1},
1903 {0, 2, -1, -1},
1904 {3, 3, 3, 3}};
1905 const int VFTetFace[][4] = {{1, 3, -1, -1},
1906 {1, 2, -1, -1},
1907 {2, 3, -1, -1},
1908 {1, 3, -1, -1},
1909 {1, 2, 1, 2}};
1910 const int VFPyrFace[][4] = {{1, 4, -1, -1},
1911 {1, 2, -1, -1},
1912 {2, 3, -1, -1},
1913 {3, 4, -1, -1},
1914 {1, 2, 3, 4}};
1915
1916 // next handle all triangular faces from tetrahedron
1917 nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
1918 for (int v = 0; v < 5; ++v)
1919 {
1920 for (int f = 0; f < nadjface[v]; ++f)
1921 {
1922 for (int i = 0; i < nfacemodes; ++i)
1923 {
1924 val = (*maxRmat[eTetrahedron])(
1925 vertMapMaxR[eTetrahedron][VFTetVert[v][f]],
1926 faceMapMaxR[eTetrahedron][VFTetFace[v][f]][i]);
1927 newmat->SetValue(vertMapMaxR[ePyramid][v],
1928 faceMapMaxR[ePyramid][VFPyrFace[v][f]][i],
1929 val);
1930 }
1931 }
1932 }
1933
1934 // Edge to face coupling
1935 // all base edges are coupled to face 0
1936 nfacemodes = (nummodesmax - 2) * (nummodesmax - 2);
1937 for (int e = 0; e < 4; ++e)
1938 {
1939 for (int i = 0; i < nummodesmax - 2; ++i)
1940 {
1941 for (int j = 0; j < nfacemodes; ++j)
1942 {
1945
1946 val = (*maxRmat[eHexahedron])(edgemapid, facemapid);
1947 newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
1948 faceMapMaxR[ePyramid][0][j], val);
1949 }
1950 }
1951 }
1952
1953 const int nadjface1[] = {1, 1, 1, 1, 2, 2, 2, 2};
1954 const int EFTetEdge[][2] = {{0, -1}, {1, -1}, {0, -1}, {2, -1},
1955 {3, 3}, {4, 4}, {5, 5}, {3, 5}};
1956 const int EFTetFace[][2] = {{1, -1}, {2, -1}, {1, -1}, {3, -1},
1957 {1, 3}, {1, 2}, {2, 3}, {1, 3}};
1958 const int EFPyrFace[][2] = {{1, -1}, {2, -1}, {3, -1}, {4, -1},
1959 {1, 4}, {1, 2}, {2, 3}, {3, 4}};
1960
1961 // next handle all triangular faces from tetrahedron
1962 nfacemodes = (nummodesmax - 3) * (nummodesmax - 2) / 2;
1963 for (int e = 0; e < 8; ++e)
1964 {
1965 for (int f = 0; f < nadjface1[e]; ++f)
1966 {
1967 for (int i = 0; i < nummodesmax - 2; ++i)
1968 {
1969 for (int j = 0; j < nfacemodes; ++j)
1970 {
1971 int edgemapid =
1972 edgeMapMaxR[eTetrahedron][EFTetEdge[e][f]][i];
1973 int facemapid =
1974 faceMapMaxR[eTetrahedron][EFTetFace[e][f]][j];
1975
1976 val = (*maxRmat[eTetrahedron])(edgemapid, facemapid);
1977 newmat->SetValue(edgeMapMaxR[ePyramid][e][i],
1978 faceMapMaxR[ePyramid][EFPyrFace[e][f]][j],
1979 val);
1980 }
1981 }
1982 }
1983 }
1984
1985 DNekScalMatSharedPtr PyrR;
1986 PyrR = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, newmat);
1987 maxRmat[ePyramid] = PyrR;
1988}
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::LibUtilities::eHexahedron, Nektar::LibUtilities::ePyramid, and Nektar::LibUtilities::eTetrahedron.
Referenced by SetUpReferenceElements().
◆ SetUpReferenceElements()
|
private |
Loop expansion and determine different variants of the transformation matrix.
Sets up multiple reference elements based on the element expansion.
Definition at line 1619 of file PreconditionerLowEnergy.cpp.
1624{
1625 std::shared_ptr<MultiRegions::ExpList> expList =
1626 ((m_linsys.lock())->GetLocMat()).lock();
1627 GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
1628
1629 LocalRegions::ExpansionSharedPtr locExp;
1630
1631 // face maps for pyramid and hybrid meshes - not need to return.
1632 map<ShapeType, Array<OneD, Array<OneD, unsigned int>>> faceMapMaxR;
1633
1634 /* Determine the maximum expansion order for all elements */
1635 int nummodesmax = 0;
1636 Array<OneD, int> Shapes(LibUtilities::SIZE_ShapeType, 0);
1637
1638 for (int n = 0; n < expList->GetNumElmts(); ++n)
1639 {
1640 locExp = expList->GetExp(n);
1641
1642 nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(0));
1643 nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(1));
1644 nummodesmax = max(nummodesmax, locExp->GetBasisNumModes(2));
1645
1646 Shapes[locExp->DetShapeType()] = 1;
1647 }
1648
1651
1653 // need Tet and Hex expansion
1654 {
1655 Shapes[eTetrahedron] = 1;
1656 Shapes[eHexahedron] = 1;
1657 }
1658
1659 StdRegions::MatrixType PreconR;
1661 {
1662 PreconR = StdRegions::ePreconRMass;
1663 }
1664 else
1665 {
1666 PreconR = StdRegions::ePreconR;
1667 }
1668
1669 Array<OneD, unsigned int> vmap;
1670 Array<OneD, Array<OneD, unsigned int>> emap;
1671 Array<OneD, Array<OneD, unsigned int>> fmap;
1672
1673 /*
1674 * Set up a transformation matrices for equal max order
1675 * polynomial meshes
1676 */
1677
1679 {
1681 // Bases for Hex element
1683 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1685 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1687 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1688
1689 // Create reference Hexahdedral expansion
1690 LocalRegions::HexExpSharedPtr HexExp;
1691
1692 HexExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(
1693 HexBa, HexBb, HexBc, hexgeom);
1694
1695 maxElmt[eHexahedron] = HexExp;
1696
1697 // Hex:
1698 HexExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1699 vertMapMaxR[eHexahedron] = vmap;
1700 edgeMapMaxR[eHexahedron] = emap;
1701 faceMapMaxR[eHexahedron] = fmap;
1702
1703 // Get hexahedral transformation matrix
1704 LocalRegions::MatrixKey HexR(PreconR, eHexahedron, *HexExp,
1705 linSysKey.GetConstFactors());
1706 maxRmat[eHexahedron] = HexExp->GetLocMatrix(HexR);
1707 }
1708
1710 {
1712 // Bases for Tetrahedral element
1714 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1716 PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
1718 PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
1719
1720 // Create reference tetrahedral expansion
1721 LocalRegions::TetExpSharedPtr TetExp;
1722
1723 TetExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(
1724 TetBa, TetBb, TetBc, tetgeom);
1725
1726 maxElmt[eTetrahedron] = TetExp;
1727
1728 TetExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1729 vertMapMaxR[eTetrahedron] = vmap;
1730 edgeMapMaxR[eTetrahedron] = emap;
1731 faceMapMaxR[eTetrahedron] = fmap;
1732
1733 // Get tetrahedral transformation matrix
1734 LocalRegions::MatrixKey TetR(PreconR, eTetrahedron, *TetExp,
1735 linSysKey.GetConstFactors());
1736 maxRmat[eTetrahedron] = TetExp->GetLocMatrix(TetR);
1737
1739 {
1740 ReSetTetMaxRMat(nummodesmax, TetExp, maxRmat, vertMapMaxR,
1741 edgeMapMaxR, faceMapMaxR);
1742 }
1743 }
1744
1746 {
1748
1749 // Bases for Pyramid element
1751 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1753 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1755 PointsKey(nummodesmax, eGaussRadauMAlpha2Beta0));
1756
1757 // Create reference pyramid expansion
1758 LocalRegions::PyrExpSharedPtr PyrExp;
1759
1760 PyrExp = MemoryManager<LocalRegions::PyrExp>::AllocateSharedPtr(
1761 PyrBa, PyrBb, PyrBc, pyrgeom);
1762
1763 maxElmt[ePyramid] = PyrExp;
1764
1765 // Pyramid:
1766 PyrExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1767 vertMapMaxR[ePyramid] = vmap;
1768 edgeMapMaxR[ePyramid] = emap;
1769 faceMapMaxR[ePyramid] = fmap;
1770
1771 // Set up Pyramid Transformation Matrix based on Tet
1772 // and Hex expansion
1773 SetUpPyrMaxRMat(nummodesmax, PyrExp, maxRmat, vertMapMaxR, edgeMapMaxR,
1774 faceMapMaxR);
1775 }
1776
1778 {
1780 // Bases for Prism element
1781 const BasisKey PrismBa(
1782 eModified_A, nummodesmax,
1783 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1784 const BasisKey PrismBb(
1785 eModified_A, nummodesmax,
1786 PointsKey(nummodesmax + 1, eGaussLobattoLegendre));
1788 PointsKey(nummodesmax, eGaussRadauMAlpha1Beta0));
1789
1790 // Create reference prismatic expansion
1791 LocalRegions::PrismExpSharedPtr PrismExp;
1792
1793 PrismExp = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(
1794 PrismBa, PrismBb, PrismBc, prismgeom);
1795 maxElmt[ePrism] = PrismExp;
1796
1797 // Prism:
1798 PrismExp->GetInverseBoundaryMaps(vmap, emap, fmap);
1799 vertMapMaxR[ePrism] = vmap;
1800 edgeMapMaxR[ePrism] = emap;
1801
1802 faceMapMaxR[ePrism] = fmap;
1803
1805 {
1806 ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
1807 edgeMapMaxR, faceMapMaxR, false);
1808 }
1809 else
1810 {
1811 // Get prismatic transformation matrix
1812 LocalRegions::MatrixKey PrismR(PreconR, ePrism, *PrismExp,
1813 linSysKey.GetConstFactors());
1814 maxRmat[ePrism] = PrismExp->GetLocMatrix(PrismR);
1815
1817 {
1818 ReSetPrismMaxRMat(nummodesmax, PrismExp, maxRmat, vertMapMaxR,
1819 edgeMapMaxR, faceMapMaxR, true);
1820 }
1821 }
1822 }
1823}
LibUtilities::CommSharedPtr m_comm
Definition: Preconditioner.h:126
SpatialDomains::HexGeomSharedPtr CreateRefHexGeom(void)
Sets up the reference hexahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1542
void ReSetPrismMaxRMat(int nummodesmax, LocalRegions::PrismExpSharedPtr &PirsmExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR, bool UseTetOnly)
Definition: PreconditionerLowEnergy.cpp:2049
void ReSetTetMaxRMat(int nummodesmax, LocalRegions::TetExpSharedPtr &TetExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
Definition: PreconditionerLowEnergy.cpp:1990
SpatialDomains::TetGeomSharedPtr CreateRefTetGeom(void)
Sets up the reference tretrahedral element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1466
SpatialDomains::PrismGeomSharedPtr CreateRefPrismGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis.
Definition: PreconditionerLowEnergy.cpp:1279
void SetUpPyrMaxRMat(int nummodesmax, LocalRegions::PyrExpSharedPtr &PyrExp, ShapeToDNekMap &maxRmat, ShapeToIntArrayMap &vertMapMaxR, ShapeToIntArrayArrayMap &edgeMapMaxR, ShapeToIntArrayArrayMap &faceMapMaxR)
Definition: PreconditionerLowEnergy.cpp:1825
SpatialDomains::PyrGeomSharedPtr CreateRefPyrGeom(void)
Sets up the reference prismatic element needed to construct a low energy basis mapping arrays.
Definition: PreconditionerLowEnergy.cpp:1376
std::shared_ptr< PrismExp > PrismExpSharedPtr
Definition: PrismExp.h:209
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::shared_ptr< HexExp > HexExpSharedPtr
Definition: HexExp.h:260
std::shared_ptr< PyrExp > PyrExpSharedPtr
Definition: PyrExp.h:185
std::shared_ptr< TetExp > TetExpSharedPtr
Definition: TetExp.h:214
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), CreateRefHexGeom(), CreateRefPrismGeom(), CreateRefPyrGeom(), CreateRefTetGeom(), Nektar::LibUtilities::eGaussLobattoLegendre, Nektar::LibUtilities::eHexahedron, Nektar::StdRegions::eMass, Nektar::LibUtilities::eModified_A, Nektar::LibUtilities::eModified_B, Nektar::LibUtilities::eModified_C, Nektar::LibUtilities::eModifiedPyr_C, Nektar::StdRegions::ePreconR, Nektar::StdRegions::ePreconRMass, Nektar::LibUtilities::ePrism, Nektar::LibUtilities::ePyramid, Nektar::LibUtilities::eTetrahedron, Nektar::MultiRegions::GlobalMatrixKey::GetConstFactors(), Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType(), Nektar::MultiRegions::Preconditioner::m_comm, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::LibUtilities::ReduceMax, ReSetPrismMaxRMat(), ReSetTetMaxRMat(), SetUpPyrMaxRMat(), and Nektar::LibUtilities::SIZE_ShapeType.
Referenced by SetupBlockTransformationMatrix().
◆ v_BuildPreconditioner()
|
overrideprotectedvirtual |
Construct the low energy preconditioner from \(\mathbf{S}_{2}\).
\[\mathbf{M}^{-1}=\left[\begin{array}{ccc} Diag[(\mathbf{S_{2}})_{vv}] & & \\ & (\mathbf{S}_{2})_{eb} & \\ & & (\mathbf{S}_{2})_{fb} \end{array}\right] \]
where \(\mathbf{R}\) is the transformation matrix and \(\mathbf{S}_{2}\) the Schur complement of the modified basis, given by
\[\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\]
where \(\mathbf{S}_{1}\) is the local schur complement matrix for each element.
- Count edges, face and add up min edges and min face sizes
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 117 of file PreconditionerLowEnergy.cpp.
118{
119 std::shared_ptr<MultiRegions::ExpList> expList =
120 ((m_linsys.lock())->GetLocMat()).lock();
121 LocalRegions::Expansion3DSharedPtr locExpansion;
122 GlobalLinSysKey linSysKey = (m_linsys.lock())->GetKey();
123
124 int i, j, k;
125 int nVerts, nEdges, nFaces;
126 int eid, fid, n, cnt, nmodes, nedgemodes, nfacemodes;
127 int nedgemodesloc;
128 NekDouble zero = 0.0;
129
130 int vMap1, vMap2, sign1, sign2;
131 int m, v, eMap1, eMap2, fMap1, fMap2;
132 int offset, globalrow, globalcol, nCoeffs;
133
134 // Periodic information
135 PeriodicMap periodicVerts;
136 PeriodicMap periodicEdges;
137 PeriodicMap periodicFaces;
138 expList->GetPeriodicEntities(periodicVerts, periodicEdges, periodicFaces);
139
140 // matrix storage
144
145 // local element static condensed matrices
146 DNekScalBlkMatSharedPtr loc_mat;
147 DNekScalMatSharedPtr bnd_mat;
148
149 DNekMatSharedPtr pRSRT;
150
151 DNekMat RS;
152 DNekMat RSRT;
153
155
156 int nDirBnd = asmMap->GetNumGlobalDirBndCoeffs();
157 int nNonDirVerts = asmMap->GetNumNonDirVertexModes();
158
159 // Vertex, edge and face preconditioner matrices
161 nNonDirVerts, nNonDirVerts, zero, vertstorage);
162
163 Array<OneD, NekDouble> vertArray(nNonDirVerts, 0.0);
164 Array<OneD, long> VertBlockToUniversalMap(nNonDirVerts, -1);
165
166 // maps for different element types
167 int n_exp = expList->GetNumElmts();
168 int nNonDirEdgeIDs = asmMap->GetNumNonDirEdges();
169 int nNonDirFaceIDs = asmMap->GetNumNonDirFaces();
170
171 set<int> edgeDirMap;
172 set<int> faceDirMap;
173 map<int, int> uniqueEdgeMap;
174 map<int, int> uniqueFaceMap;
175
176 // this should be of size total number of local edges + faces
177 Array<OneD, int> modeoffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs, 0);
178 Array<OneD, int> globaloffset(1 + nNonDirEdgeIDs + nNonDirFaceIDs, 0);
179
180 const Array<OneD, const ExpListSharedPtr> &bndCondExp =
181 expList->GetBndCondExpansions();
182 LocalRegions::Expansion2DSharedPtr bndCondFaceExp;
183 const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
184 &bndConditions = expList->GetBndConditions();
185
186 int meshVertId;
187 int meshEdgeId;
188 int meshFaceId;
189
190 const Array<OneD, const int> &extradiredges = asmMap->GetExtraDirEdges();
191 for (i = 0; i < extradiredges.size(); ++i)
192 {
193 meshEdgeId = extradiredges[i];
194 edgeDirMap.insert(meshEdgeId);
195 }
196
197 // Determine which boundary edges and faces have dirichlet values
198 for (i = 0; i < bndCondExp.size(); i++)
199 {
200 cnt = 0;
201 for (j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
202 {
203 bndCondFaceExp =
204 std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
205 bndCondExp[i]->GetExp(j));
206 if (bndConditions[i]->GetBoundaryConditionType() ==
208 {
209 for (k = 0; k < bndCondFaceExp->GetNtraces(); k++)
210 {
211 meshEdgeId = bndCondFaceExp->GetGeom()->GetEid(k);
212 if (edgeDirMap.count(meshEdgeId) == 0)
213 {
214 edgeDirMap.insert(meshEdgeId);
215 }
216 }
217 meshFaceId = bndCondFaceExp->GetGeom()->GetGlobalID();
218 faceDirMap.insert(meshFaceId);
219 }
220 }
221 }
222
223 int dof = 0;
224 int maxFaceDof = 0;
225 int maxEdgeDof = 0;
226 int nlocalNonDirEdges = 0;
227 int nlocalNonDirFaces = 0;
228 int ntotalentries = 0;
229
230 map<int, int> EdgeSize;
231 map<int, int> FaceSize;
232 map<int, pair<int, int>> FaceModes;
233
234 /// - Count edges, face and add up min edges and min face sizes
235 for (n = 0; n < n_exp; ++n)
236 {
237 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
238
239 nEdges = locExpansion->GetNedges();
240 for (j = 0; j < nEdges; ++j)
241 {
242 int nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
243 meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
244 ->GetGeom3D()
245 ->GetEid(j);
246 if (EdgeSize.count(meshEdgeId) == 0)
247 {
248 EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
249 }
250 else
251 {
252 EdgeSize[meshEdgeId] =
253 min(EdgeSize[meshEdgeId], nEdgeInteriorCoeffs);
254 }
255 }
256
257 nFaces = locExpansion->GetNtraces();
258 for (j = 0; j < nFaces; ++j)
259 {
260 int nFaceInteriorCoeffs = locExpansion->GetTraceIntNcoeffs(j);
261 meshFaceId = locExpansion->GetGeom3D()->GetFid(j);
262
263 if (FaceSize.count(meshFaceId) == 0)
264 {
265 FaceSize[meshFaceId] = nFaceInteriorCoeffs;
266
267 int m0, m1;
268 locExpansion->GetTraceNumModes(j, m0, m1,
269 locExpansion->GetTraceOrient(j));
270 FaceModes[meshFaceId] = pair<int, int>(m0, m1);
271 }
272 else
273 {
274 if (nFaceInteriorCoeffs < FaceSize[meshFaceId])
275 {
276 FaceSize[meshFaceId] = nFaceInteriorCoeffs;
277 int m0, m1;
278 locExpansion->GetTraceNumModes(
279 j, m0, m1, locExpansion->GetTraceOrient(j));
280 FaceModes[meshFaceId] = pair<int, int>(m0, m1);
281 }
282 }
283 }
284 }
285
286 bool verbose = expList->GetSession()->DefinesCmdLineArgument("verbose");
287
288 // For parallel runs need to check have minimum of edges and faces over
289 // partition boundaries
291 {
292 int EdgeSizeLen = EdgeSize.size();
293 int FaceSizeLen = FaceSize.size();
294 Array<OneD, long> FacetMap(EdgeSizeLen + FaceSizeLen, -1);
295 Array<OneD, NekDouble> FacetLen(EdgeSizeLen + FaceSizeLen, -1);
296
297 map<int, int>::iterator it;
298
299 cnt = 0;
300 int maxid = 0;
301 for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
302 {
303 FacetMap[cnt] = it->first;
304 maxid = max(it->first, maxid);
305 FacetLen[cnt] = it->second;
306 }
307 maxid++;
308
310
311 for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
312 {
313 FacetMap[cnt] = it->first + maxid;
314 FacetLen[cnt] = it->second;
315 }
316
317 // Exchange vertex data over different processes
320
321 cnt = 0;
322 for (it = EdgeSize.begin(); it != EdgeSize.end(); ++it, ++cnt)
323 {
324 it->second = (int)FacetLen[cnt];
325 }
326
327 for (it = FaceSize.begin(); it != FaceSize.end(); ++it, ++cnt)
328 {
329 it->second = (int)FacetLen[cnt];
330 }
331 }
332
333 // Loop over all the elements in the domain and compute total edge
334 // DOF and set up unique ordering.
335 map<int, int> nblks;
336 int matrixlocation = 0;
337
338 // First do periodic edges
339 for (auto &pIt : periodicEdges)
340 {
341 meshEdgeId = pIt.first;
342
343 if (edgeDirMap.count(meshEdgeId) == 0)
344 {
345 dof = EdgeSize[meshEdgeId];
346 if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
347 {
348 bool SetUpNewEdge = true;
349
350 for (i = 0; i < pIt.second.size(); ++i)
351 {
352 if (!pIt.second[i].isLocal)
353 {
354 continue;
355 }
356
357 int meshEdgeId2 = pIt.second[i].id;
358 if (edgeDirMap.count(meshEdgeId2) == 0)
359 {
360 if (uniqueEdgeMap.count(meshEdgeId2) != 0)
361 {
362 // set unique map to same location
363 uniqueEdgeMap[meshEdgeId] =
364 uniqueEdgeMap[meshEdgeId2];
365 SetUpNewEdge = false;
366 }
367 }
368 else
369 {
370 edgeDirMap.insert(meshEdgeId);
371 SetUpNewEdge = false;
372 }
373 }
374
375 if (SetUpNewEdge)
376 {
377 uniqueEdgeMap[meshEdgeId] = matrixlocation;
378 globaloffset[matrixlocation] += ntotalentries;
379 modeoffset[matrixlocation] = dof * dof;
380 ntotalentries += dof * dof;
381 nblks[matrixlocation++] = dof;
382 }
383 }
384 }
385 }
386
387 for (cnt = n = 0; n < n_exp; ++n)
388 {
389 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
390
391 for (j = 0; j < locExpansion->GetNedges(); ++j)
392 {
393 meshEdgeId = locExpansion->GetGeom()->GetEid(j);
394 dof = EdgeSize[meshEdgeId];
395 maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
396
397 if (edgeDirMap.count(meshEdgeId) == 0)
398 {
399 if (uniqueEdgeMap.count(meshEdgeId) == 0 && dof > 0)
400
401 {
402 uniqueEdgeMap[meshEdgeId] = matrixlocation;
403
404 globaloffset[matrixlocation] += ntotalentries;
405 modeoffset[matrixlocation] = dof * dof;
406 ntotalentries += dof * dof;
407 nblks[matrixlocation++] = dof;
408 }
409 nlocalNonDirEdges += dof * dof;
410 }
411 }
412 }
413
414 // Loop over all the elements in the domain and compute max face
415 // DOF. Reduce across all processes to get universal maximum.
416 // - Periodic faces
417 for (auto &pIt : periodicFaces)
418 {
419 meshFaceId = pIt.first;
420
421 if (faceDirMap.count(meshFaceId) == 0)
422 {
423 dof = FaceSize[meshFaceId];
424
425 if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
426 {
427 bool SetUpNewFace = true;
428
429 if (pIt.second[0].isLocal)
430 {
431 int meshFaceId2 = pIt.second[0].id;
432
433 if (faceDirMap.count(meshFaceId2) == 0)
434 {
435 if (uniqueFaceMap.count(meshFaceId2) != 0)
436 {
437 // set unique map to same location
438 uniqueFaceMap[meshFaceId] =
439 uniqueFaceMap[meshFaceId2];
440 SetUpNewFace = false;
441 }
442 }
443 else // set face to be a Dirichlet face
444 {
445 faceDirMap.insert(meshFaceId);
446 SetUpNewFace = false;
447 }
448 }
449
450 if (SetUpNewFace)
451 {
452 uniqueFaceMap[meshFaceId] = matrixlocation;
453
454 modeoffset[matrixlocation] = dof * dof;
455 globaloffset[matrixlocation] += ntotalentries;
456 ntotalentries += dof * dof;
457 nblks[matrixlocation++] = dof;
458 }
459 }
460 }
461 }
462
463 for (cnt = n = 0; n < n_exp; ++n)
464 {
465 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
466
467 for (j = 0; j < locExpansion->GetNtraces(); ++j)
468 {
469 meshFaceId = locExpansion->GetGeom()->GetFid(j);
470
471 dof = FaceSize[meshFaceId];
472 maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
473
474 if (faceDirMap.count(meshFaceId) == 0)
475 {
476 if (uniqueFaceMap.count(meshFaceId) == 0 && dof > 0)
477 {
478 uniqueFaceMap[meshFaceId] = matrixlocation;
479 modeoffset[matrixlocation] = dof * dof;
480 globaloffset[matrixlocation] += ntotalentries;
481 ntotalentries += dof * dof;
482 nblks[matrixlocation++] = dof;
483 }
484 nlocalNonDirFaces += dof * dof;
485 }
486 }
487 }
488
491
492 // Allocate arrays for block to universal map (number of expansions * p^2)
493 Array<OneD, long> BlockToUniversalMap(ntotalentries, -1);
494 Array<OneD, int> localToGlobalMatrixMap(
495 nlocalNonDirEdges + nlocalNonDirFaces, -1);
496
497 // Allocate arrays to store matrices (number of expansions * p^2)
498 Array<OneD, NekDouble> BlockArray(nlocalNonDirEdges + nlocalNonDirFaces,
499 0.0);
500
501 int matrixoffset = 0;
502 int vGlobal;
503 int uniEdgeOffset = 0;
504
505 // Need to obtain a fixed offset for the universal number
506 // of the faces which come after the edge numbering
507 for (n = 0; n < n_exp; ++n)
508 {
509 for (j = 0; j < locExpansion->GetNedges(); ++j)
510 {
511 meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
512 ->GetGeom3D()
513 ->GetEid(j);
514
515 uniEdgeOffset = max(meshEdgeId, uniEdgeOffset);
516 }
517 }
518 uniEdgeOffset++;
519
521 uniEdgeOffset = uniEdgeOffset * maxEdgeDof * maxEdgeDof;
522
523 for (n = 0; n < n_exp; ++n)
524 {
525 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
526
527 // loop over the edges of the expansion
528 for (j = 0; j < locExpansion->GetNedges(); ++j)
529 {
530 // get mesh edge id
531 meshEdgeId = locExpansion->GetGeom()->GetEid(j);
532
533 nedgemodes = EdgeSize[meshEdgeId];
534
535 if (edgeDirMap.count(meshEdgeId) == 0 && nedgemodes > 0)
536 {
537 // Determine the Global edge offset
538 int edgeOffset = globaloffset[uniqueEdgeMap[meshEdgeId]];
539
540 // Determine a universal map offset
541 int uniOffset = meshEdgeId;
542 auto pIt = periodicEdges.find(meshEdgeId);
543 if (pIt != periodicEdges.end())
544 {
545 for (int l = 0; l < pIt->second.size(); ++l)
546 {
547 uniOffset = min(uniOffset, pIt->second[l].id);
548 }
549 }
550 uniOffset = uniOffset * maxEdgeDof * maxEdgeDof;
551
552 for (k = 0; k < nedgemodes * nedgemodes; ++k)
553 {
554 vGlobal = edgeOffset + k;
555 localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
556 BlockToUniversalMap[vGlobal] = uniOffset + k + 1;
557 }
558 matrixoffset += nedgemodes * nedgemodes;
559 }
560 }
561
562 Array<OneD, unsigned int> faceInteriorMap;
563 Array<OneD, int> faceInteriorSign;
564 // loop over the faces of the expansion
565 for (j = 0; j < locExpansion->GetNtraces(); ++j)
566 {
567 // get mesh face id
568 meshFaceId = locExpansion->GetGeom()->GetFid(j);
569
570 nfacemodes = FaceSize[meshFaceId];
571
572 // Check if face has dirichlet values
573 if (faceDirMap.count(meshFaceId) == 0 && nfacemodes > 0)
574 {
575 // Determine the Global edge offset
576 int faceOffset = globaloffset[uniqueFaceMap[meshFaceId]];
577 // Determine a universal map offset
578 int uniOffset = meshFaceId;
579 // use minimum face edge when periodic
580 auto pIt = periodicFaces.find(meshFaceId);
581 if (pIt != periodicFaces.end())
582 {
583 uniOffset = min(uniOffset, pIt->second[0].id);
584 }
585 uniOffset = uniOffset * maxFaceDof * maxFaceDof;
586
587 for (k = 0; k < nfacemodes * nfacemodes; ++k)
588 {
589 vGlobal = faceOffset + k;
590
591 localToGlobalMatrixMap[matrixoffset + k] = vGlobal;
592
593 BlockToUniversalMap[vGlobal] =
594 uniOffset + uniEdgeOffset + k + 1;
595 }
596 matrixoffset += nfacemodes * nfacemodes;
597 }
598 }
599 }
600
601 matrixoffset = 0;
602
603 map<int, int>::iterator it;
604 Array<OneD, unsigned int> n_blks(nblks.size() + 1);
605 n_blks[0] = nNonDirVerts;
606 for (i = 1, it = nblks.begin(); it != nblks.end(); ++it)
607 {
608 n_blks[i++] = it->second;
609 }
610
612 blkmatStorage);
613
614 // Here we loop over the expansion and build the block low energy
615 // preconditioner as well as the block versions of the transformation
616 // matrices.
617 for (cnt = n = 0; n < n_exp; ++n)
618 {
619 locExpansion = expList->GetExp(n)->as<LocalRegions::Expansion3D>();
620 nCoeffs = locExpansion->NumBndryCoeffs();
621
622 // Get correct transformation matrix for element type
624
625 pRSRT = MemoryManager<DNekMat>::AllocateSharedPtr(nCoeffs, nCoeffs,
626 zero, storage);
627 RSRT = (*pRSRT);
628
629 nVerts = locExpansion->GetGeom()->GetNumVerts();
630 nEdges = locExpansion->GetGeom()->GetNumEdges();
631 nFaces = locExpansion->GetGeom()->GetNumFaces();
632
633 // Get statically condensed matrix
634 loc_mat = (m_linsys.lock())->GetStaticCondBlock(n);
635
636 // Extract boundary block (elemental S1)
637 bnd_mat = loc_mat->GetBlock(0, 0);
638
639 // offset by number of rows
640 offset = bnd_mat->GetRows();
641
642 DNekScalMat &S = (*bnd_mat);
643
644 DNekMat Sloc(nCoeffs, nCoeffs);
645
646 // For variable p we need to just use the active modes
647 NekDouble val;
648
649 for (int i = 0; i < nCoeffs; ++i)
650 {
651 for (int j = 0; j < nCoeffs; ++j)
652 {
653 val = S(i, j) * m_variablePmask[cnt + i] *
654 m_variablePmask[cnt + j];
655 Sloc.SetValue(i, j, val);
656 }
657 }
658
659 // Calculate R*S*trans(R)
660 RSRT = R * Sloc * Transpose(R);
661
662 // loop over vertices of the element and return the vertex map
663 // for each vertex
664 for (v = 0; v < nVerts; ++v)
665 {
666 vMap1 = locExpansion->GetVertexMap(v);
667
668 // Get vertex map
669 globalrow = asmMap->GetLocalToGlobalBndMap(cnt + vMap1) - nDirBnd;
670
671 if (globalrow >= 0)
672 {
673 for (m = 0; m < nVerts; ++m)
674 {
675 vMap2 = locExpansion->GetVertexMap(m);
676
677 // global matrix location (without offset due to
678 // dirichlet values)
679 globalcol =
680 asmMap->GetLocalToGlobalBndMap(cnt + vMap2) - nDirBnd;
681
682 // offset for dirichlet conditions
683 if (globalcol == globalrow)
684 {
685 // modal connectivity between elements
686 sign1 = asmMap->GetLocalToGlobalBndSign(cnt + vMap1);
687 sign2 = asmMap->GetLocalToGlobalBndSign(cnt + vMap2);
688
689 vertArray[globalrow] +=
690 sign1 * sign2 * RSRT(vMap1, vMap2);
691
692 meshVertId = locExpansion->GetGeom()->GetVid(v);
693
694 auto pIt = periodicVerts.find(meshVertId);
695 if (pIt != periodicVerts.end())
696 {
697 for (k = 0; k < pIt->second.size(); ++k)
698 {
699 meshVertId = min(meshVertId, pIt->second[k].id);
700 }
701 }
702
703 VertBlockToUniversalMap[globalrow] = meshVertId + 1;
704 }
705 }
706 }
707 }
708
709 // loop over edges of the element and return the edge map
710 for (eid = 0; eid < nEdges; ++eid)
711 {
712 meshEdgeId = locExpansion->as<LocalRegions::Expansion3D>()
713 ->GetGeom3D()
714 ->GetEid(eid);
715
716 nedgemodes = EdgeSize[meshEdgeId];
717 if (nedgemodes)
718 {
719 nedgemodesloc = locExpansion->GetEdgeNcoeffs(eid) - 2;
720 DNekMatSharedPtr m_locMat =
722 nedgemodes, nedgemodes, zero, storage);
723
724 Array<OneD, unsigned int> edgemodearray =
725 locExpansion->GetEdgeInverseBoundaryMap(eid);
726
727 if (edgeDirMap.count(meshEdgeId) == 0)
728 {
729 for (v = 0; v < nedgemodesloc; ++v)
730 {
731 eMap1 = edgemodearray[v];
732 sign1 = asmMap->GetLocalToGlobalBndSign(cnt + eMap1);
733
734 if (sign1 == 0)
735 {
736 continue;
737 }
738
739 for (m = 0; m < nedgemodesloc; ++m)
740 {
741 eMap2 = edgemodearray[m];
742
743 // modal connectivity between elements
744 sign2 =
745 asmMap->GetLocalToGlobalBndSign(cnt + eMap2);
746
747 NekDouble globalEdgeValue =
748 sign1 * sign2 * RSRT(eMap1, eMap2);
749
750 if (sign2 != 0)
751 {
752 // if(eMap1 == eMap2)
753 BlockArray[matrixoffset + v * nedgemodes + m] =
754 globalEdgeValue;
755 }
756 }
757 }
758 matrixoffset += nedgemodes * nedgemodes;
759 }
760 }
761 }
762
763 // loop over faces of the element and return the face map
764 for (fid = 0; fid < nFaces; ++fid)
765 {
766 meshFaceId = locExpansion->as<LocalRegions::Expansion3D>()
767 ->GetGeom3D()
768 ->GetFid(fid);
769
770 nfacemodes = FaceSize[meshFaceId];
771 if (nfacemodes > 0)
772 {
773 DNekMatSharedPtr m_locMat =
775 nfacemodes, nfacemodes, zero, storage);
776
777 if (faceDirMap.count(meshFaceId) == 0)
778 {
779 Array<OneD, unsigned int> facemodearray;
780 StdRegions::Orientation faceOrient =
781 locExpansion->GetTraceOrient(fid);
782
783 auto pIt = periodicFaces.find(meshFaceId);
784 if (pIt != periodicFaces.end())
785 {
786 if (meshFaceId == min(meshFaceId, pIt->second[0].id))
787 {
788 faceOrient = DeterminePeriodicFaceOrient(
789 faceOrient, pIt->second[0].orient);
790 }
791 }
792
793 facemodearray = locExpansion->GetTraceInverseBoundaryMap(
794 fid, faceOrient, FaceModes[meshFaceId].first,
795 FaceModes[meshFaceId].second);
796
797 for (v = 0; v < nfacemodes; ++v)
798 {
799 fMap1 = facemodearray[v];
800
801 sign1 = asmMap->GetLocalToGlobalBndSign(cnt + fMap1);
802
803 ASSERTL1(sign1 != 0,
804 "Something is wrong since we "
805 "shoudl only be extracting modes for "
806 "lowest order expansion");
807
808 for (m = 0; m < nfacemodes; ++m)
809 {
810 fMap2 = facemodearray[m];
811
812 // modal connectivity between elements
813 sign2 =
814 asmMap->GetLocalToGlobalBndSign(cnt + fMap2);
815
816 ASSERTL1(sign2 != 0,
817 "Something is wrong since "
818 "we shoudl only be extracting modes for "
819 "lowest order expansion");
820
821 // Get the face-face value from the
822 // low energy matrix (S2)
823 NekDouble globalFaceValue =
824 sign1 * sign2 * RSRT(fMap1, fMap2);
825
826 // local face value to global face value
827 // if(fMap1 == fMap2)
828 BlockArray[matrixoffset + v * nfacemodes + m] =
829 globalFaceValue;
830 }
831 }
832 matrixoffset += nfacemodes * nfacemodes;
833 }
834 }
835 }
836
837 // offset for the expansion
838 cnt += nCoeffs;
839 }
840
841 // Exchange vertex data over different processes
844
845 Array<OneD, NekDouble> GlobalBlock(ntotalentries, 0.0);
846 if (ntotalentries)
847 {
848 // Assemble edge matrices of each process
849 Vmath::Assmb(BlockArray.size(), BlockArray, localToGlobalMatrixMap,
850 GlobalBlock);
851 }
852
853 // Exchange edge & face data over different processes
856
857 // Populate vertex block
858 for (int i = 0; i < nNonDirVerts; ++i)
859 {
860 VertBlk->SetValue(i, i, 1.0 / vertArray[i]);
861 }
862
863 // Set the first block to be the diagonal of the vertex space
864 m_BlkMat->SetBlock(0, 0, VertBlk);
865
866 // Build the edge and face matrices from the vector
867 DNekMatSharedPtr gmat;
868
869 offset = 0;
870 // -1 since we ignore vert block
872 {
873 nmodes = n_blks[1 + loc];
874 gmat = MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
875 storage);
876
877 for (v = 0; v < nmodes; ++v)
878 {
879 for (m = 0; m < nmodes; ++m)
880 {
881 NekDouble Value = GlobalBlock[offset + v * nmodes + m];
882 gmat->SetValue(v, m, Value);
883 }
884 }
886 offset += modeoffset[loc];
887 }
888
889 // invert blocks.
891 for (i = 1; i < totblks; ++i)
892 {
894 if (nmodes)
895 {
896 DNekMatSharedPtr tmp_mat =
897 MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, nmodes, zero,
898 storage);
899
900 tmp_mat = m_BlkMat->GetBlock(i, i);
901 tmp_mat->Invert();
902
903 m_BlkMat->SetBlock(i, i, tmp_mat);
904 }
905 }
906}
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
DNekBlkMatSharedPtr m_BlkMat
Definition: PreconditionerLowEnergy.h:81
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:270
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:192
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:48
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Definition: MultiRegions.hpp:177
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
Definition: AssemblyMapCG.cpp:2271
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:67
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
Definition: BlockMatrix.hpp:235
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:857
Definition: GsLib.hpp:151
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL1, Vmath::Assmb(), Nektar::MultiRegions::DeterminePeriodicFaceOrient(), Nektar::eDIAGONAL, Nektar::SpatialDomains::eDirichlet, Nektar::eFULL, Gs::Gather(), Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs(), Nektar::LocalRegions::Expansion::GetGeom(), Nektar::StdRegions::StdExpansion3D::GetNedges(), Gs::gs_add, Gs::gs_min, Gs::Init(), CG_Iterations::loc, m_BlkMat, Nektar::MultiRegions::Preconditioner::m_comm, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_variablePmask, Nektar::StdRegions::StdExpansion::NumBndryCoeffs(), Nektar::LibUtilities::ReduceMax, and Nektar::Transpose().
◆ v_DoPreconditioner()
|
overrideprotectedvirtual |
Apply the low energy preconditioner during the conjugate gradient routine
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 912 of file PreconditionerLowEnergy.cpp.
916{
918 "set up for Global iteratives sovles");
921 int nNonDir = nGlobal - nDir;
922 DNekBlkMat &M = (*m_BlkMat);
923
924 NekVector<NekDouble> r(nNonDir, pInput, eWrapper);
925
927
928 z = M * r;
929}
std::vector< double > z(NPUPPER)
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, BlockMatrixTag > DNekBlkMat
Definition: NekTypeDefs.hpp:58
References ASSERTL0, Nektar::eWrapper, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and Nektar::UnitTests::z().
◆ v_DoTransformBasisFromLowEnergy()
|
overrideprotectedvirtual |
Multiply by the block inverse transformation matrix This transforms the bassi from Low Energy to original basis.
Note; not typically required
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 1127 of file PreconditionerLowEnergy.cpp.
1129{
1131
1132 ASSERTL1(pInput.size() >= nLocBndDofs,
1133 "Input array is smaller than nLocBndDofs");
1134 ASSERTL1(pOutput.size() >= nLocBndDofs,
1135 "Output array is smaller than nLocBndDofs");
1136
1137 // Block inverse transformation matrix
1139
1140 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
1141
1142 // Multiply by the inverse transformation matrix
1143 int cnt = 0;
1144 int cnt1 = 0;
1146 {
1150 &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1151 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pOutput.get() + cnt,
1152 nbndcoeffs);
1153 cnt += nbndcoeffs * nexp;
1154 cnt1 += nexp;
1155 }
1156}
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:385
References ASSERTL1, Blas::Dgemm(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and m_sameBlock.
◆ v_DoTransformBasisToLowEnergy()
|
overrideprotectedvirtual |
Transform the basis vector to low energy space.
As the conjugate gradient system is solved for the low energy basis, the solution vector \(\mathbf{x}\) must be transformed to the low energy basis i.e. \(\overline{\mathbf{x}}=\mathbf{R}\mathbf{x}\).
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 1051 of file PreconditionerLowEnergy.cpp.
1053{
1055
1056 // Block transformation matrix
1058
1059 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
1060
1061 // Apply mask in case of variable P
1063
1064 // Multiply by the block transformation matrix
1065 int cnt = 0;
1066 int cnt1 = 0;
1068 {
1072 &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1073 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pInOut.get() + cnt,
1074 nbndcoeffs);
1075 cnt += nbndcoeffs * nexp;
1076 cnt1 += nexp;
1077 }
1078}
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:207
References Blas::Dgemm(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, m_variablePmask, and Vmath::Vmul().
◆ v_DoTransformCoeffsFromLowEnergy()
|
overrideprotectedvirtual |
transform the solution coeffiicents from low energy back to the original coefficient space.
After the conjugate gradient routine the output vector is in the low energy basis and must be trasnformed back to the original basis in order to get the correct solution out. the solution vector i.e. \(\mathbf{x}=\mathbf{R^{T}}\mathbf{\overline{x}}\).
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 1089 of file PreconditionerLowEnergy.cpp.
1091{
1093
1094 ASSERTL1(pInOut.size() >= nLocBndDofs,
1095 "Output array is not greater than the nLocBndDofs");
1096
1097 // Block transposed transformation matrix
1099
1100 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInOut.get());
1101
1102 // Multiply by the transpose of block transformation matrix
1103 int cnt = 0;
1104 int cnt1 = 0;
1106 {
1110 &(R.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1111 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pInOut.get() + cnt,
1112 nbndcoeffs);
1113 cnt += nbndcoeffs * nexp;
1114 cnt1 += nexp;
1115 }
1116
1118}
References ASSERTL1, Blas::Dgemm(), Nektar::MultiRegions::Preconditioner::m_locToGloMap, m_RBlk, m_sameBlock, m_variablePmask, and Vmath::Vmul().
◆ v_DoTransformCoeffsToLowEnergy()
|
overrideprotectedvirtual |
Multiply by the block tranposed inverse transformation matrix (R^T)^{-1} which is equivlaent to transforming coefficients to LowEnergy space.
In JCP 2001 paper on low energy this is seen as (C^T)^{-1}
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 1165 of file PreconditionerLowEnergy.cpp.
1167{
1169
1170 ASSERTL1(pInput.size() >= nLocBndDofs,
1171 "Input array is less than nLocBndDofs");
1172 ASSERTL1(pOutput.size() >= nLocBndDofs,
1173 "Output array is less than nLocBndDofs");
1174
1175 // Block inverse transformation matrix
1177
1178 Array<OneD, NekDouble> pLocalIn(nLocBndDofs, pInput.get());
1179
1180 // Multiply by the transpose of block transformation matrix
1181 int cnt = 0;
1182 int cnt1 = 0;
1184 {
1188 &(invR.GetBlock(cnt1, cnt1)->GetPtr()[0]), nbndcoeffs,
1189 pLocalIn.get() + cnt, nbndcoeffs, 0.0, pOutput.get() + cnt,
1190 nbndcoeffs);
1191 cnt += nbndcoeffs * nexp;
1192 cnt1 += nexp;
1193 }
1194}
References ASSERTL1, Blas::Dgemm(), m_InvRBlk, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and m_sameBlock.
◆ v_InitObject()
|
overrideprotectedvirtual |
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 72 of file PreconditionerLowEnergy.cpp.
73{
75
77 solvertype == ePETScStaticCond,
78 "Solver type not valid");
79
80 std::shared_ptr<MultiRegions::ExpList> expList =
81 ((m_linsys.lock())->GetLocMat()).lock();
82
83 m_comm = expList->GetComm();
84
85 LocalRegions::ExpansionSharedPtr locExpansion;
86
87 locExpansion = expList->GetExp(0);
88
89 int nDim = locExpansion->GetShapeDimension();
90
92
93 // Set up block transformation matrix
95
96 // Sets up variable p mask
97 CreateVariablePMask();
98}
void SetupBlockTransformationMatrix(void)
Definition: PreconditionerLowEnergy.cpp:936
void CreateVariablePMask(void)
Definition: PreconditionerLowEnergy.cpp:1252
References ASSERTL0, CreateVariablePMask(), Nektar::MultiRegions::eIterativeStaticCond, Nektar::MultiRegions::ePETScStaticCond, Nektar::MultiRegions::Preconditioner::m_comm, Nektar::MultiRegions::Preconditioner::m_linsys, Nektar::MultiRegions::Preconditioner::m_locToGloMap, and SetupBlockTransformationMatrix().
◆ v_TransformedSchurCompl()
|
overrideprotectedvirtual |
Set up the transformed block matrix system.
Sets up a block elemental matrix in which each of the block matrix is the low energy equivalent i.e. \(\mathbf{S}_{2}=\mathbf{R}\mathbf{S}_{1}\mathbf{R}^{T}\)
Reimplemented from Nektar::MultiRegions::Preconditioner.
Definition at line 1203 of file PreconditionerLowEnergy.cpp.
1205{
1206 std::shared_ptr<MultiRegions::ExpList> expList =
1207 ((m_linsys.lock())->GetLocMat()).lock();
1208
1209 LocalRegions::ExpansionSharedPtr locExpansion;
1210 locExpansion = expList->GetExp(n);
1211 unsigned int nbnd = locExpansion->NumBndryCoeffs();
1212
1214 DNekMatSharedPtr pS2 =
1215 MemoryManager<DNekMat>::AllocateSharedPtr(nbnd, nbnd, 0.0, storage);
1216
1217 // transformation matrices
1219
1220 // Original Schur Complement
1221 DNekScalMat &S1 = (*loc_mat);
1222
1223 DNekMat Sloc(nbnd, nbnd);
1224
1225 // For variable p we need to just use the active modes
1226 NekDouble val;
1227
1228 for (int i = 0; i < nbnd; ++i)
1229 {
1230 for (int j = 0; j < nbnd; ++j)
1231 {
1232 val = S1(i, j) * m_variablePmask[bndoffset + i] *
1233 m_variablePmask[bndoffset + j];
1234 Sloc.SetValue(i, j, val);
1235 }
1236 }
1237
1238 // create low energy matrix
1239 DNekMat &S2 = (*pS2);
1240
1241 S2 = R * Sloc * Transpose(R);
1242
1243 DNekScalMatSharedPtr return_val;
1244 return_val = MemoryManager<DNekScalMat>::AllocateSharedPtr(1.0, pS2);
1245
1246 return return_val;
1247}
References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::eFULL, Nektar::MultiRegions::Preconditioner::m_linsys, m_RBlk, m_variablePmask, and Nektar::Transpose().
Member Data Documentation
◆ className
|
static |
Initial value:
=
"LowEnergy Preconditioning")
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
static PreconditionerSharedPtr create(const std::shared_ptr< GlobalLinSys > &plinsys, const std::shared_ptr< AssemblyMap > &pLocToGloMap)
Creates an instance of this class.
Definition: PreconditionerLowEnergy.h:57
PreconFactory & GetPreconFactory()
Definition: Preconditioner.cpp:68
Name of class.
Registers the class with the Factory.
Definition at line 69 of file PreconditionerLowEnergy.h.
◆ m_BlkMat
|
protected |
Definition at line 81 of file PreconditionerLowEnergy.h.
Referenced by v_BuildPreconditioner().
◆ m_InvRBlk
|
protected |
Definition at line 83 of file PreconditionerLowEnergy.h.
Referenced by SetupBlockTransformationMatrix(), v_DoTransformBasisFromLowEnergy(), and v_DoTransformCoeffsToLowEnergy().
◆ m_RBlk
|
protected |
Definition at line 82 of file PreconditionerLowEnergy.h.
Referenced by SetupBlockTransformationMatrix(), v_BuildPreconditioner(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_TransformedSchurCompl().
◆ m_sameBlock
|
protected |
Definition at line 91 of file PreconditionerLowEnergy.h.
Referenced by SetupBlockTransformationMatrix(), v_DoTransformBasisFromLowEnergy(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_DoTransformCoeffsToLowEnergy().
◆ m_signChange
|
protected |
Definition at line 87 of file PreconditionerLowEnergy.h.
Referenced by CreateVariablePMask().
◆ m_variablePmask
Definition at line 85 of file PreconditionerLowEnergy.h.
Referenced by CreateVariablePMask(), v_BuildPreconditioner(), v_DoTransformBasisToLowEnergy(), v_DoTransformCoeffsFromLowEnergy(), and v_TransformedSchurCompl().
