Commit 567a7102 authored by Andreas Wagner's avatar Andreas Wagner
Browse files

starts simple tutorial for heat equation

parent 9c2d8ec7
Pipeline #32143 failed with stages
in 112 minutes and 49 seconds
waLBerla_add_executable( NAME HeatEquation
FILES HeatEquation.cpp
DEPENDS hyteg core)
waLBerla_link_files_to_builddir( *.prm )
/*
* Copyright (c) 2017-2019 Dominik Bartuschat, Dominik Thoennes, Marcus Mohr, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <cmath>
#include <core/math/Constants.h>
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/mpi/MPIManager.h"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/gridtransferoperators/P1toP1LinearProlongation.hpp"
#include "hyteg/gridtransferoperators/P1toP1LinearRestriction.hpp"
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/solvers/CGSolver.hpp"
#include "hyteg/solvers/GeometricMultigridSolver.hpp"
#include "hyteg/solvers/SymmetricGaussSeidelSmoother.hpp"
using walberla::real_t;
using walberla::uint_c;
using walberla::uint_t;
template < typename OpType >
std::shared_ptr< hyteg::Solver< OpType > > create_solver( const std::shared_ptr< hyteg::PrimitiveStorage >& storage,
const walberla::Config::BlockHandle& parameters )
{
const uint_t minLevel = parameters.getParameter< uint_t >( "minLevel" );
const uint_t maxLevel = parameters.getParameter< uint_t >( "maxLevel" );
const uint_t max_coarse_iter = parameters.getParameter< uint_t >( "max_coarse_iter" );
const real_t coarse_tolerance = parameters.getParameter< real_t >( "coarse_tolerance" );
const real_t outer_tolerance = parameters.getParameter< real_t >( "outer_tolerance" );
auto smoother = std::make_shared< hyteg::SymmetricGaussSeidelSmoother< OpType > >();
auto coarseGridSolver =
std::make_shared< hyteg::CGSolver< OpType > >( storage, minLevel, minLevel, max_coarse_iter, coarse_tolerance );
auto restrictionOperator = std::make_shared< hyteg::P1toP1LinearRestriction >();
auto prolongationOperator = std::make_shared< hyteg::P1toP1LinearProlongation >();
auto multiGridSolver = std::make_shared< hyteg::GeometricMultigridSolver< OpType > >(
storage, smoother, coarseGridSolver, restrictionOperator, prolongationOperator, minLevel, maxLevel );
auto outerSolver =
std::make_shared< hyteg::CGSolver< OpType > >( storage, minLevel, maxLevel, 1000, outer_tolerance, multiGridSolver );
outerSolver->setPrintInfo( true );
return outerSolver;
}
hyteg::P1LinearCombinationForm create_lhs_form( double dt )
{
using MassFormType = hyteg::P1FenicsForm< p1_mass_cell_integral_0_otherwise, p1_tet_mass_cell_integral_0_otherwise >;
using LaplaceFormType =
hyteg::P1FenicsForm< p1_diffusion_cell_integral_0_otherwise, p1_tet_diffusion_cell_integral_0_otherwise >;
hyteg::P1LinearCombinationForm lhsForm;
lhsForm.addOwnedForm< MassFormType >( 1. );
lhsForm.addOwnedForm< LaplaceFormType >( dt );
return lhsForm;
}
real_t analytic_solution( double t, const hyteg::Point3D& x )
{
constexpr auto pi = walberla::math::pi;
const double k_x = 1.;
const double k_y = 2;
return std::exp( -pi * pi * ( k_x * k_x + k_y * k_y ) * t ) * std::cos( pi * k_x * x[0] ) * std::cos( pi * k_y * x[1] );
}
std::shared_ptr< hyteg::PrimitiveStorage > create_storage()
{
hyteg::Point2D lowerRight( { -1, -1 } );
hyteg::Point2D upperLeft( { +1, +1 } );
auto meshInfo = hyteg::MeshInfo::meshRectangle( lowerRight, upperLeft, hyteg::MeshInfo::CROSS, 2, 2 );
auto setupStorage = std::make_shared< hyteg::SetupPrimitiveStorage >(
meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
// Neumann BC everywhere
setupStorage->setMeshBoundaryFlagsOnBoundary( 2, 0, true );
return std::make_shared< hyteg::PrimitiveStorage >( *setupStorage );
}
int main( int argc, char** argv )
{
walberla::Environment env( argc, argv );
walberla::mpi::MPIManager::instance()->useWorldComm();
walberla::shared_ptr< walberla::config::Config > cfg( new walberla::config::Config );
cfg->readParameterFile( "./HeatEquation.prm" );
walberla::Config::BlockHandle parameters = cfg->getOneBlock( "Parameters" );
parameters.listParameters();
const uint_t minLevel = parameters.getParameter< uint_t >( "minLevel" );
const uint_t maxLevel = parameters.getParameter< uint_t >( "maxLevel" );
auto storage = create_storage();
hyteg::P1Function< real_t > residual( "residual", storage, minLevel, maxLevel );
hyteg::P1Function< real_t > error( "error", storage, minLevel, maxLevel );
hyteg::P1Function< real_t > rightHandSide( "rightHandSide", storage, minLevel, maxLevel );
hyteg::P1Function< real_t > u_now( "u_now", storage, minLevel, maxLevel );
hyteg::P1Function< real_t > u_analytic( "u_analytic", storage, minLevel, maxLevel );
hyteg::P1Function< real_t > laplaceTimesFunction( "laplaceTimesFunction", storage, minLevel, maxLevel );
double t = 0;
double t_end = 1e-1;
auto analytic_solution_currently = [&t]( const hyteg::Point3D& p ) { return analytic_solution( t, p ); };
u_now.interpolate( analytic_solution_currently, maxLevel, hyteg::All );
double dt = 1e-3;
hyteg::P1ConstantLinearCombinationOperator lhsOp( storage, minLevel, maxLevel, create_lhs_form( dt ) );
hyteg::P1ConstantMassOperator rhsOp( storage, minLevel, maxLevel );
auto solver = create_solver< hyteg::P1ConstantLinearCombinationOperator >( storage, parameters );
hyteg::VTKOutput vtkOutput( ".", "HeatEquation", storage );
if ( parameters.getParameter< bool >( "vtkOutput" ) )
{
vtkOutput.add( u_now );
vtkOutput.add( residual );
vtkOutput.add( rightHandSide );
vtkOutput.add( u_analytic );
}
const auto num_time_steps = static_cast< uint_t >( std::ceil( t_end / dt ) );
for ( uint_t t_idx = 0; t_idx < num_time_steps; t_idx += 1 )
{
t += dt;
rhsOp.apply( u_now, rightHandSide, maxLevel, hyteg::All, hyteg::Replace );
solver->solve( lhsOp, u_now, rightHandSide, maxLevel );
{
u_analytic.interpolate( analytic_solution_currently, maxLevel, hyteg::All );
error.assign( { 1.0, -1.0 }, { u_now, u_analytic }, maxLevel, hyteg::Inner );
real_t errorNorm = std::sqrt( error.dotGlobal( error, maxLevel, hyteg::Inner ) );
WALBERLA_LOG_INFO_ON_ROOT( "error = " << errorNorm )
}
if ( parameters.getParameter< bool >( "vtkOutput" ) )
{
vtkOutput.write( maxLevel, t_idx );
}
}
}
Parameters
{
minLevel 2;
maxLevel 6;
max_outer_iter 20;
max_coarse_iter 10000;
mg_tolerance 1e-16;
outer_tolerance 1e-12;
coarse_tolerance 1e-16;
vtkOutput true;
}
......@@ -17,3 +17,4 @@ waLBerla_add_executable( NAME 04_Indexing
add_subdirectory(05_FullAppP1GMG)
add_subdirectory(06_FullAppPlumeInCube)
add_subdirectory(07_IsoviscousConvectionAnnulus)
add_subdirectory(08_HeatEquation)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment