Skip to content

GitLab

Commit 474976d2 authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Merge branch 'kohl/dg-3D-rebased-01' into 'master' 

DG implementation

See merge request !502
parents 47f72935 7ebe6dc1
Pipeline #38972 passed with stages
in 136 minutes and 19 seconds
Expand all Hide whitespace changes
Inline Side-by-side
apps/Chorin.cpp
View file @ 474976d2
......@@ -19,7 +19,7 @@
*/
#include <core/math/Constants.h>
#include "hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include "hyteg/dgfunctionspace_old/DGUpwindOperator.hpp"
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
......@@ -121,11 +121,11 @@ int main( int argc, char* argv[] )
hyteg::P1Function< real_t > res( "res", storage, minLevel, maxLevel );
hyteg::P1Function< real_t > ones( "ones", storage, minLevel, maxLevel );
auto u_dg = std::make_shared< hyteg::DGFunction< real_t > >( "u_dg", storage, minLevel, maxLevel );
auto v_dg = std::make_shared< hyteg::DGFunction< real_t > >( "v_dg", storage, minLevel, maxLevel );
auto u_dg = std::make_shared< hyteg::DGFunction_old< real_t > >( "u_dg", storage, minLevel, maxLevel );
auto v_dg = std::make_shared< hyteg::DGFunction_old< real_t > >( "v_dg", storage, minLevel, maxLevel );
auto u_dg_old = std::make_shared< hyteg::DGFunction< real_t > >( "u_dg", storage, minLevel, maxLevel );
auto v_dg_old = std::make_shared< hyteg::DGFunction< real_t > >( "v_dg", storage, minLevel, maxLevel );
auto u_dg_old = std::make_shared< hyteg::DGFunction_old< real_t > >( "u_dg", storage, minLevel, maxLevel );
auto v_dg_old = std::make_shared< hyteg::DGFunction_old< real_t > >( "v_dg", storage, minLevel, maxLevel );
hyteg::P1ConstantLaplaceOperator A( storage, minLevel, maxLevel );
hyteg::P1ConstantLaplaceOperator Ascaled( storage, minLevel, maxLevel );
......
apps/Surrogates/VariableCoeffLaplace_Surrogates.cpp
View file @ 474976d2
......@@ -384,11 +384,11 @@ void showStencilFunction(std::shared_ptr<PrimitiveStorage> storage, const uint_t
P2Form_divKgrad form;
form.setGeometryMap(face.getGeometryMap());
Point3D x0(face.coords[0]), x;
Point3D x0(face.getCoordinates()[0]), x;
real_t h = 1.0 / (walberla::real_c(levelinfo::num_microvertices_per_edge(level) - 1));
Point3D d0 = h * (face.coords[1] - face.coords[0]);
Point3D d2 = h * (face.coords[2] - face.coords[0]);
Point3D d0 = h * (face.getCoordinates()[1] - face.getCoordinates()[0]);
Point3D d2 = h * (face.getCoordinates()[2] - face.getCoordinates()[0]);
// directions
const Point3D dirS = -d2;
......@@ -412,8 +412,8 @@ void showStencilFunction(std::shared_ptr<PrimitiveStorage> storage, const uint_t
P2::variablestencil::macroface::assembleStencil(form, x, dirS, dirSE, dirE, dirN, dirNW, dirW, dirNE,
VtVStencil, EtVStencil, VtEStencil, EtEStencil);
uint_t i = it.col();
uint_t j = it.row();
idx_t i = it.col();
idx_t j = it.row();
// VERTEX DoF
if (!vertexdof::macroface::isVertexOnBoundary(level, it))
......
apps/benchmarks/ApplyPerformanceAnalysis-2D-P2/AbstractApply.cpp
View file @ 474976d2
......@@ -37,9 +37,9 @@ void apply_2d_vertex_to_vertex( double* WALBERLA_RESTRICT dstPtr,
const uint_t level )
{
uint_t rowsize = levelinfo::num_microvertices_per_edge( level );
for ( uint_t j = 1; j < rowsize - 2; ++j )
for ( idx_t j = 1; j < idx_t( rowsize ) - 2; ++j )
{
for ( uint_t i = 1; i < rowsize -j - 1; ++i )
for ( idx_t i = 1; i < idx_t( rowsize ) -j - 1; ++i )
{
auto tmp = real_t( 0 );
for ( const auto direction : vertexdof::macroface::neighborsWithCenter )
......
apps/benchmarks/CompareInterpolate.cpp
View file @ 474976d2
......@@ -46,17 +46,17 @@ inline void interpolateStdFunction( Face&
uint_t rowsize = levelinfo::num_microvertices_per_edge( Level );
Point3D x, x0;
auto dstPtr = faceMemory->getPointer( Level );
x0 = face.coords[0];
Point3D d0 = ( face.coords[1] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.coords[2] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
x0 = face.getCoordinates()[0];
Point3D d0 = ( face.getCoordinates()[1] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.getCoordinates()[2] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
uint_t inner_rowsize = rowsize;
for( uint_t i = 1; i < rowsize - 2; ++i )
for( idx_t i = 1; i < idx_t( rowsize ) - 2; ++i )
{
x = x0;
x += real_c( i ) * d2 + d0;
for( uint_t j = 1; j < inner_rowsize - 2; ++j )
for( idx_t j = 1; j < idx_t (inner_rowsize ) - 2; ++j )
{
dstPtr[vertexdof::macroface::indexFromVertex( Level, j, i, stencilDirection::VERTEX_C )] = expr( x );
x += d0;
......@@ -74,17 +74,17 @@ inline void
uint_t rowsize = levelinfo::num_microvertices_per_edge( Level );
Point3D x, x0;
auto dstPtr = faceMemory->getPointer( Level );
x0 = face.coords[0];
Point3D d0 = ( face.coords[1] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.coords[2] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
x0 = face.getCoordinates()[0];
Point3D d0 = ( face.getCoordinates()[1] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.getCoordinates()[2] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
uint_t inner_rowsize = rowsize;
for( uint_t i = 1; i < rowsize - 2; ++i )
for( idx_t i = 1; i < idx_t( rowsize ) - 2; ++i )
{
x = x0;
x += real_c( i ) * d2 + d0;
for( uint_t j = 1; j < inner_rowsize - 2; ++j )
for( idx_t j = 1; j < idx_t( inner_rowsize ) - 2; ++j )
{
dstPtr[vertexdof::macroface::indexFromVertex( Level, j, i, stencilDirection::VERTEX_C )] = expr( x );
x += d0;
......@@ -103,17 +103,17 @@ inline void interpolateFunctor( Face&
uint_t rowsize = levelinfo::num_microvertices_per_edge( Level );
Point3D x, x0;
auto dstPtr = faceMemory->getPointer( Level );
x0 = face.coords[0];
Point3D d0 = ( face.coords[1] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.coords[2] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
x0 = face.getCoordinates()[0];
Point3D d0 = ( face.getCoordinates()[1] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.getCoordinates()[2] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
uint_t inner_rowsize = rowsize;
for( uint_t i = 1; i < rowsize - 2; ++i )
for( idx_t i = 1; i < idx_t( rowsize ) - 2; ++i )
{
x = x0;
x += real_c( i ) * d2 + d0;
for( uint_t j = 1; j < inner_rowsize - 2; ++j )
for( idx_t j = 1; j < idx_t( inner_rowsize ) - 2; ++j )
{
dstPtr[vertexdof::macroface::indexFromVertex( Level, j, i, stencilDirection::VERTEX_C )] = exprFunctor( x );
x += d0;
......@@ -130,18 +130,18 @@ inline void interpolateWithoutFunction( Face& face, const PrimitiveDataID< Funct
uint_t rowsize = levelinfo::num_microvertices_per_edge( Level );
Point3D x, x0;
auto dstPtr = faceMemory->getPointer( Level );
x0 = face.coords[0];
Point3D d0 = ( face.coords[1] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.coords[2] - face.coords[0] ) / ( walberla::real_c( rowsize - 1 ) );
x0 = face.getCoordinates()[0];
Point3D d0 = ( face.getCoordinates()[1] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
Point3D d2 = ( face.getCoordinates()[2] - face.getCoordinates()[0] ) / ( walberla::real_c( rowsize - 1 ) );
uint_t inner_rowsize = rowsize;
for( uint_t i = 1; i < rowsize - 2; ++i )
for( idx_t i = 1; i < idx_t( rowsize ) - 2; ++i )
{
x = x0;
x += real_c( i ) * d2 + d0;
for( uint_t j = 1; j < inner_rowsize - 2; ++j )
for( idx_t j = 1; j < idx_t( inner_rowsize ) - 2; ++j )
{
dstPtr[vertexdof::macroface::indexFromVertex( Level, j, i, stencilDirection::VERTEX_C )] =
sqrt( x[0] * x[0] + x[1] * x[1] );
......
apps/dg0_transport.cpp
View file @ 474976d2
......@@ -21,8 +21,8 @@
#include <core/timing/Timer.h>
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/dgfunctionspace_old/DGUpwindOperator.hpp"
#include "hyteg/dgfunctionspace_old/DGFunction.hpp"
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
......@@ -74,8 +74,8 @@ int main(int argc, char* argv[])
std::shared_ptr< hyteg::PrimitiveStorage> storage = std::make_shared< hyteg::PrimitiveStorage>(setupStorage);
std::shared_ptr< hyteg::DGFunction<real_t>> c_old = std::make_shared< hyteg::DGFunction<real_t>>("c_old", storage, minLevel, maxLevel);
std::shared_ptr< hyteg::DGFunction<real_t>> c = std::make_shared< hyteg::DGFunction<real_t>>("c", storage, minLevel, maxLevel);
std::shared_ptr< hyteg::DGFunction_old<real_t>> c_old = std::make_shared< hyteg::DGFunction_old<real_t>>("c_old", storage, minLevel, maxLevel);
std::shared_ptr< hyteg::DGFunction_old<real_t>> c = std::make_shared< hyteg::DGFunction_old<real_t>>("c", storage, minLevel, maxLevel);
std::shared_ptr< hyteg::P1Function<real_t>> u = std::make_shared< hyteg::P1Function<real_t>>("u", storage, minLevel, maxLevel);
std::shared_ptr< hyteg::P1Function<real_t>> v = std::make_shared< hyteg::P1Function<real_t>>("v", storage, minLevel, maxLevel);
......
apps/geophysics_plume.cpp
View file @ 474976d2
......@@ -24,8 +24,8 @@
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/P1P1StokesOperator.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/dgfunctionspace/DGFunction.hpp"
#include "hyteg/dgfunctionspace/DGUpwindOperator.hpp"
#include "hyteg/dgfunctionspace_old/DGFunction.hpp"
#include "hyteg/dgfunctionspace_old/DGUpwindOperator.hpp"
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/gridtransferoperators/P1P1StokesToP1P1StokesProlongation.hpp"
#include "hyteg/gridtransferoperators/P1P1StokesToP1P1StokesRestriction.hpp"
......@@ -101,10 +101,10 @@ int main( int argc, char* argv[] )
#endif
// Setting up Functions
auto c_old = std::make_shared< hyteg::DGFunction< real_t > >( "c", storage, minLevel, maxLevel );
auto c = std::make_shared< hyteg::DGFunction< real_t > >( "c", storage, minLevel, maxLevel );
auto c_old = std::make_shared< hyteg::DGFunction_old< real_t > >( "c", storage, minLevel, maxLevel );
auto c = std::make_shared< hyteg::DGFunction_old< real_t > >( "c", storage, minLevel, maxLevel );
auto f_dg = std::make_shared< hyteg::DGFunction< real_t > >( "f_dg", storage, minLevel, maxLevel );
auto f_dg = std::make_shared< hyteg::DGFunction_old< real_t > >( "f_dg", storage, minLevel, maxLevel );
auto r = std::make_shared< hyteg::P1StokesFunction< real_t > >( "r", storage, minLevel, maxLevel );
auto f = std::make_shared< hyteg::P1StokesFunction< real_t > >( "f", storage, minLevel, maxLevel );
......@@ -165,11 +165,11 @@ int main( int argc, char* argv[] )
for ( uint_t lvl = minLevel; lvl <= maxLevel; ++lvl )
{
uint_t tmpDofStokes = numberOfGlobalDoFs< hyteg::P1StokesFunctionTag >( *storage, lvl );
uint_t tmpDofTemp = numberOfGlobalDoFs< hyteg::DGFunctionTag >( *storage, lvl );
// uint_t tmpDofTemp = numberOfGlobalDoFs< hyteg::DGFunctionTag >( *storage, lvl );
WALBERLA_LOG_DETAIL_ON_ROOT( "Stokes DoFs on level " << lvl << " : " << tmpDofStokes );
WALBERLA_LOG_DETAIL_ON_ROOT( "Temperature DoFs on level " << lvl << " : " << tmpDofTemp );
// WALBERLA_LOG_DETAIL_ON_ROOT( "Temperature DoFs on level " << lvl << " : " << tmpDofTemp );
totalGlobalDofsStokes += tmpDofStokes;
totalGlobalDofsTemp += tmpDofTemp;
// totalGlobalDofsTemp += tmpDofTemp;
}
WALBERLA_LOG_INFO_ON_ROOT( "Total Stokes DoFs on all level :" << totalGlobalDofsStokes );
WALBERLA_LOG_INFO_ON_ROOT( "Total Temperature DoFs on all level :" << totalGlobalDofsTemp );
......
src/coupling_hyteg_convection_particles/MMOCTransport.hpp
View file @ 474976d2
......@@ -58,39 +58,40 @@ enum class TimeSteppingScheme
RK4, // fourth order
};
const static std::vector< std::vector< real_t > > RK_A_ExplicitEuler = {{{0}}};
const static std::vector< real_t > RK_b_ExplicitEuler = {1};
const static std::vector< real_t > RK_c_ExplicitEuler = {0};
const static std::vector< std::vector< real_t > > RK_A_ExplicitEuler = { { { 0 } } };
const static std::vector< real_t > RK_b_ExplicitEuler = { 1 };
const static std::vector< real_t > RK_c_ExplicitEuler = { 0 };
const static std::vector< std::vector< real_t > > RK_A_RK3 = {{{0, 0, 0}, {0.5, 0, 0}, {-1., 2., 0}}};
const static std::vector< real_t > RK_b_RK3 = {1. / 6., 2. / 3., 1. / 6.};
const static std::vector< real_t > RK_c_RK3 = {0, 0.5, 1.};
const static std::vector< std::vector< real_t > > RK_A_RK3 = { { { 0, 0, 0 }, { 0.5, 0, 0 }, { -1., 2., 0 } } };
const static std::vector< real_t > RK_b_RK3 = { 1. / 6., 2. / 3., 1. / 6. };
const static std::vector< real_t > RK_c_RK3 = { 0, 0.5, 1. };
const static std::vector< std::vector< real_t > > RK_A_Ralston = {{{0, 0, 0}, {0.5, 0, 0}, {0, 0.75, 0}}};
const static std::vector< real_t > RK_b_Ralston = {2. / 9., 1. / 3., 4. / 9.};
const static std::vector< real_t > RK_c_Ralston = {0, 0.5, 0.75};
const static std::vector< std::vector< real_t > > RK_A_Ralston = { { { 0, 0, 0 }, { 0.5, 0, 0 }, { 0, 0.75, 0 } } };
const static std::vector< real_t > RK_b_Ralston = { 2. / 9., 1. / 3., 4. / 9. };
const static std::vector< real_t > RK_c_Ralston = { 0, 0.5, 0.75 };
const static std::vector< std::vector< real_t > > RK_A_RK4 = {{{0, 0, 0, 0}, {0.5, 0, 0, 0}, {0, 0.5, 0, 0}, {0, 0, 1, 0}}};
const static std::vector< real_t > RK_b_RK4 = {1. / 6., 1. / 3., 1. / 3., 1. / 6.};
const static std::vector< real_t > RK_c_RK4 = {0, 0.5, 0.5, 1.};
const static std::vector< std::vector< real_t > > RK_A_RK4 = {
{ { 0, 0, 0, 0 }, { 0.5, 0, 0, 0 }, { 0, 0.5, 0, 0 }, { 0, 0, 1, 0 } } };
const static std::vector< real_t > RK_b_RK4 = { 1. / 6., 1. / 3., 1. / 3., 1. / 6. };
const static std::vector< real_t > RK_c_RK4 = { 0, 0.5, 0.5, 1. };
const static std::map< TimeSteppingScheme, std::vector< std::vector< real_t > > > RK_A = {
{TimeSteppingScheme::ExplicitEuler, RK_A_ExplicitEuler},
{TimeSteppingScheme::RK3, RK_A_RK3},
{TimeSteppingScheme::Ralston, RK_A_Ralston},
{TimeSteppingScheme::RK4, RK_A_RK4}};
{ TimeSteppingScheme::ExplicitEuler, RK_A_ExplicitEuler },
{ TimeSteppingScheme::RK3, RK_A_RK3 },
{ TimeSteppingScheme::Ralston, RK_A_Ralston },
{ TimeSteppingScheme::RK4, RK_A_RK4 } };
const static std::map< TimeSteppingScheme, std::vector< real_t > > RK_b = {
{TimeSteppingScheme::ExplicitEuler, RK_b_ExplicitEuler},
{TimeSteppingScheme::RK3, RK_b_RK3},
{TimeSteppingScheme::Ralston, RK_b_Ralston},
{TimeSteppingScheme::RK4, RK_b_RK4}};
{ TimeSteppingScheme::ExplicitEuler, RK_b_ExplicitEuler },
{ TimeSteppingScheme::RK3, RK_b_RK3 },
{ TimeSteppingScheme::Ralston, RK_b_Ralston },
{ TimeSteppingScheme::RK4, RK_b_RK4 } };
const static std::map< TimeSteppingScheme, std::vector< real_t > > RK_c = {
{TimeSteppingScheme::ExplicitEuler, RK_c_ExplicitEuler},
{TimeSteppingScheme::RK3, RK_c_RK3},
{TimeSteppingScheme::Ralston, RK_c_Ralston},
{TimeSteppingScheme::RK4, RK_c_RK4}};
{ TimeSteppingScheme::ExplicitEuler, RK_c_ExplicitEuler },
{ TimeSteppingScheme::RK3, RK_c_RK3 },
{ TimeSteppingScheme::Ralston, RK_c_Ralston },
{ TimeSteppingScheme::RK4, RK_c_RK4 } };
inline std::vector< PrimitiveID > getNeighboringPrimitives( const PrimitiveID& primitiveID, const PrimitiveStorage& storage )
{
......@@ -147,12 +148,12 @@ inline void updateParticlePosition( const PrimitiveStorage&
Point3D computationalLocation;
face->getGeometryMap()->evalFinv( toPoint3D( p->getPosition() ), computationalLocation );
Point2D computationalLocation2D( {computationalLocation[0], computationalLocation[1]} );
Point2D computationalLocation2D( { computationalLocation[0], computationalLocation[1] } );
if ( isPointInTriangle( computationalLocation2D,
Point2D( {face->getCoordinates().at( 0 )[0], face->getCoordinates().at( 0 )[1]} ),
Point2D( {face->getCoordinates().at( 1 )[0], face->getCoordinates().at( 1 )[1]} ),
Point2D( {face->getCoordinates().at( 2 )[0], face->getCoordinates().at( 2 )[1]} ) ) )
Point2D( { face->getCoordinates().at( 0 )[0], face->getCoordinates().at( 0 )[1] } ),
Point2D( { face->getCoordinates().at( 1 )[0], face->getCoordinates().at( 1 )[1] } ),
Point2D( { face->getCoordinates().at( 2 )[0], face->getCoordinates().at( 2 )[1] } ) ) )
{
p->setContainingPrimitive( faceID );
foundByPointLocation = true;
......@@ -168,13 +169,13 @@ inline void updateParticlePosition( const PrimitiveStorage&
const auto neighborFace = storage.getFace( neighborFaceID );
Point3D computationalLocationNeighbor;
neighborFace->getGeometryMap()->evalFinv( toPoint3D( p->getPosition() ), computationalLocationNeighbor );
Point2D computationalLocationNeighbor2D( {computationalLocationNeighbor[0], computationalLocationNeighbor[1]} );
Point2D computationalLocationNeighbor2D( { computationalLocationNeighbor[0], computationalLocationNeighbor[1] } );
if ( isPointInTriangle(
computationalLocationNeighbor2D,
Point2D( {neighborFace->getCoordinates().at( 0 )[0], neighborFace->getCoordinates().at( 0 )[1]} ),
Point2D( {neighborFace->getCoordinates().at( 1 )[0], neighborFace->getCoordinates().at( 1 )[1]} ),
Point2D( {neighborFace->getCoordinates().at( 2 )[0], neighborFace->getCoordinates().at( 2 )[1]} ) ) )
Point2D( { neighborFace->getCoordinates().at( 0 )[0], neighborFace->getCoordinates().at( 0 )[1] } ),
Point2D( { neighborFace->getCoordinates().at( 1 )[0], neighborFace->getCoordinates().at( 1 )[1] } ),
Point2D( { neighborFace->getCoordinates().at( 2 )[0], neighborFace->getCoordinates().at( 2 )[1] } ) ) )
{
// set it to the first neighbor we found to contain the particle
p->setContainingPrimitive( neighborFaceID );
......@@ -262,8 +263,7 @@ inline void updateParticlePosition( const PrimitiveStorage&
else
{
// check for neighbor cells if we did not find it in its previous cell
const auto& neighboringCells = cell->getIndirectNeighborCellIDs();
for ( const auto& neighborCellID : neighboringCells )
for ( const auto& neighborCellID : cell->getIndirectNeighborCellIDsOverVertices() )
{
WALBERLA_ASSERT( storage.cellExistsLocally( neighborCellID ) ||
storage.cellExistsInNeighborhood( neighborCellID ) );
......@@ -308,8 +308,7 @@ inline void updateParticlePosition( const PrimitiveStorage&
}
else
{
const auto& neighboringCells = cell->getIndirectNeighborCellIDs();
for ( const auto& neighborCellID : neighboringCells )
for ( const auto& neighborCellID : cell->getIndirectNeighborCellIDsOverVertices() )
{
WALBERLA_ASSERT( storage.cellExistsLocally( neighborCellID ) ||
storage.cellExistsInNeighborhood( neighborCellID ) );
......@@ -524,10 +523,10 @@ inline uint_t initializeParticles( walberla::convection_particles::data::Particl
particleStorage.clear();
const uint_t rank = uint_c( walberla::mpi::MPIManager::instance()->rank() );
const uint_t rank = uint_c( walberla::mpi::MPIManager::instance()->rank() );
// const std::vector< std::vector< real_t > >& A = RK_A.at( timeSteppingScheme ); //this seems to be unused?
const std::vector< real_t >& b = RK_b.at( timeSteppingScheme );
const uint_t rkStages = b.size();
const std::vector< real_t >& b = RK_b.at( timeSteppingScheme );
const uint_t rkStages = b.size();
if constexpr ( std::is_same< FunctionType, P1Function< real_t > >::value )
{
......@@ -806,10 +805,10 @@ inline void particleIntegration( walberla::convection_particles::data::ParticleS
storage.getTimingTree()->start( "Evaluate at particle position" );
std::vector< real_t > results( {0, 0} );
std::vector< real_t > resultsLastTimeStep( {0, 0} );
std::vector< FunctionType > functions = {ux, uy};
std::vector< FunctionType > functionsLastTimeStep = {uxLastTimeStep, uyLastTimeStep};
std::vector< real_t > results( { 0, 0 } );
std::vector< real_t > resultsLastTimeStep( { 0, 0 } );
std::vector< FunctionType > functions = { ux, uy };
std::vector< FunctionType > functionsLastTimeStep = { uxLastTimeStep, uyLastTimeStep };
if ( storage.hasGlobalCells() )
{
results.push_back( 0 );
......@@ -944,15 +943,22 @@ inline void evaluateTemperature( walberla::convection_particles::data::ParticleS
const int TAG = 98234;
#ifdef _MSC_VER
//need a first receive to avoid blocking on windows while communicating with itself
int selfCommMessage = 0;
int selfCommMessage = 0;
MPI_Request selfCommRequest;
MPI_Irecv(&selfCommMessage, 1, MPI_INT, walberla::mpi::MPIManager::instance()->rank(), TAG, walberla::mpi::MPIManager::instance()->comm(), &selfCommRequest);
MPI_Irecv( &selfCommMessage,
1,
MPI_INT,
walberla::mpi::MPIManager::instance()->rank(),
TAG,
walberla::mpi::MPIManager::instance()->comm(),
&selfCommRequest );
#endif
for ( const auto& p : particleStorage )
{
if ( numParticlesToSendToRank.count( p.getStartProcess() ) == 0 ){
if ( numParticlesToSendToRank.count( p.getStartProcess() ) == 0 )
{
numParticlesToSendToRank[p.getStartProcess()] = 0;
sendRequests[p.getStartProcess()] = MPI_Request();
sendRequests[p.getStartProcess()] = MPI_Request();
}
numParticlesToSendToRank[p.getStartProcess()]++;
}
......@@ -977,21 +983,23 @@ inline void evaluateTemperature( walberla::convection_particles::data::ParticleS
}
int numReceivedParticleLocations = 0;
#ifdef _MSC_VER
#ifdef _MSC_VER
//get self communication
{
int ready;
int ready;
MPI_Status status;
MPI_Test(&selfCommRequest, &ready, &status);
if(ready){
numReceivedParticleLocations = selfCommMessage;
MPI_Test( &selfCommRequest, &ready, &status );
if ( ready )
{
numReceivedParticleLocations = selfCommMessage;
numParticlesToReceiveFromRank[uint_c( status.MPI_SOURCE )] = numReceivedParticleLocations;
}else{
WALBERLA_LOG_INFO("somethings very wrong here");
}
else
{
WALBERLA_LOG_INFO( "somethings very wrong here" );
}
}
#endif
#endif
while ( numReceivedParticleLocations < numberOfCreatedParticles )
{
MPI_Status status;
......@@ -1184,16 +1192,16 @@ class MMOCTransport
particleLocationRadius_ = 0.1 * MeshQuality::getMinimalEdgeLength( storage, maxLevel );
}
void step( const FunctionType& c,
const vecfun_t& u,
const vecfun_t& uLastTimeStep,
const uint_t& level,
const DoFType& flag,
const real_t& dt,
const uint_t& innerSteps,
const bool& resetParticles = true,
const bool& globalMaxLimiter = true,
const bool& setParticlesOutsideDomainToZero = false )
void step( const FunctionType& c,
const vecfun_t& u,
const vecfun_t& uLastTimeStep,
const uint_t& level,
const DoFType& flag,
const real_t& dt,
const uint_t& innerSteps,
const bool& resetParticles = true,
const bool& globalMaxLimiter = true,
const bool& setParticlesOutsideDomainToZero = false )
{
uint_t aux = u.getDimension() == 3 ? 2 : 0;
step( c,
......@@ -1213,18 +1221,18 @@ class MMOCTransport
}
template < typename MassOperator >
void step( const FunctionType& c,
const vecfun_t& u,
const vecfun_t& uLastTimeStep,
const uint_t& level,
const DoFType& flag,
const real_t& dt,
const uint_t& innerSteps,
const MassOperator& massOperator,
const real_t& allowedRelativeMassDifference,
const real_t& dtPertubationAdjustedAdvection,
const bool& globalMaxLimiter = true,
const bool& setParticlesOutsideDomainToZero = false )
void step( const FunctionType& c,
const vecfun_t& u,
const vecfun_t& uLastTimeStep,
const uint_t& level,
const DoFType& flag,
const real_t& dt,
const uint_t& innerSteps,
const MassOperator& massOperator,
const real_t& allowedRelativeMassDifference,
const real_t& dtPertubationAdjustedAdvection,
const bool& globalMaxLimiter = true,
const bool& setParticlesOutsideDomainToZero = false )
{
uint_t aux = u.getDimension() == 3 ? 2 : 0;
step( c,
......@@ -1271,7 +1279,7 @@ class MMOCTransport
if ( resetParticles )
{
cOld_.assign( {1.0}, {c}, level, All );
cOld_.assign( { 1.0 }, { c }, level, All );
storage_->getTimingTree()->start( "Particle initialization" );
numberOfCreatedParticles_ =
initializeParticles( particleStorage_, *storage_, c, ux, uy, uz, level, Inner, timeSteppingSchemeConvection_, 0 );
......@@ -1335,8 +1343,8 @@ class MMOCTransport
cMinus_.copyBoundaryConditionFromFunction( c );
cAdjusted_.copyBoundaryConditionFromFunction( c );
cPlus_.assign( {1.0}, {c}, level, All );
cMinus_.assign( {1.0}, {c}, level, All );
cPlus_.assign( { 1.0 }, { c }, level, All );
cMinus_.assign( { 1.0 }, { c }, level, All );
// calculate old mass
massOperator.apply( c, cTmp_, level, flag );
......@@ -1407,11 +1415,11 @@ class MMOCTransport
if ( massAfter <= massBefore )
{
cAdjusted_.interpolate( maxAssignment, {cPlus_, cMinus_}, level );
cAdjusted_.interpolate( maxAssignment, { cPlus_, cMinus_ }, level );
}
else
{
cAdjusted_.interpolate( minAssignment, {cPlus_, cMinus_}, level );
cAdjusted_.interpolate( minAssignment, { cPlus_, cMinus_ }, level );
}
// calculate adjustment mass
......@@ -1420,7 +1428,7 @@ class MMOCTransport
auto theta = ( massBefore - massAdjusted ) / ( massAfter - massAdjusted );
c.assign( {theta, 1 - theta}, {c, cAdjusted_}, level, flag );
c.assign( { theta, 1 - theta }, { c, cAdjusted_ }, level, flag );
}
const std::shared_ptr< PrimitiveStorage > storage_;
......
src/hyteg/CMakeLists.txt
View file @ 474976d2
......@@ -35,6 +35,8 @@ add_subdirectory( boundary )
add_subdirectory( composites )
add_subdirectory( solvers )
add_subdirectory( dgfunctionspace )
add_subdirectory( dgfunctionspace_old )
add_subdirectory( volumedofspace )
add_subdirectory( mixedoperators )
add_subdirectory( eigen )
add_subdirectory( elementwiseoperators )
......@@ -50,6 +52,7 @@ add_subdirectory( petsc )
add_subdirectory( geometry )
add_subdirectory( indexing )
add_subdirectory( communication )
add_subdirectory( p0functionspace )
add_subdirectory( p1functionspace )
add_subdirectory( facedofspace )
add_subdirectory( facedofspace_old )
add_subdirectory( polynomial )
src/hyteg/Levelinfo.hpp
View file @ 474976d2
/*
* Copyright (c) 2017-2019 Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-2022 Dominik Thoennes, Nils Kohl.