Transport equation EOC is working. (a8b0cd34) · Commits · TNL / tnl-dev

examples/transport-equation/tnl-transport-equation-eoc.h

+11 −1

Original line number	Diff line number	Diff line
		@@ -41,19 +41,29 @@ template< typename ConfigTag >class advectionConfig
		static void configSetup( Config::ConfigDescription& config )
		{
		config.addDelimiter( "Transport equation settings:" );
		config.addDelimiter( "Initial condition" );
		config.addEntry< String >( "initial-condition", "Set type of initial condition.", "heaviside-sphere" );
		config.addEntryEnum< String >( "heaviside-sphere" );
		config.addEntry < double >( "x-center", "x-center for paraboloids.", 0.0 );
		config.addEntry < double >( "y-center", "y-center for paraboloids.", 0.0 );
		config.addEntry < double >( "z-center", "z-center for paraboloids.", 0.0 );
		config.addEntry < double >( "coefficient", "Coefficient for paraboloids.", -1.0 );
		config.addEntry < double >( "radius", "Radius for paraboloids.", 1.0 );

		config.addDelimiter( "Velocity field" );
		config.addEntry< String >( "velocity-field", "Type of velocity field.", "constant" );
		config.addEntryEnum< String >( "constant" );
		Functions::VectorField< 3, Functions::Analytic::Constant< 3 > >::configSetup( config, "velocity-field-" );

		config.addDelimiter( "Numerical scheme" );
		typedef Meshes::Grid< 3 > MeshType;
		Operators::Advection::LaxFridrichs< MeshType >::configSetup( config );

		config.addDelimiter( "Boundary conditions" );
		config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
		config.addEntryEnum< String >( "dirichlet" );
		config.addEntryEnum< String >( "neumann" );
		config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
		config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions.", 0.0 );
		}
		};

examples/transport-equation/transportEquationProblemEoc.h

+5 −2

Original line number	Diff line number	Diff line
		@@ -32,7 +32,7 @@ public transportEquationProblem< Mesh, BoundaryCondition, RightHandSide, Differe
		typedef typename Mesh::DeviceType DeviceType;
		typedef typename DifferentialOperator::IndexType IndexType;
		typedef Functions::MeshFunction< Mesh > MeshFunctionType;
		typedef PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
		typedef transportEquationProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator > BaseType;
		typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
		typedef SharedPointer< DifferentialOperator > DifferentialOperatorPointer;
		typedef SharedPointer< BoundaryCondition > BoundaryConditionPointer;
		@@ -40,6 +40,7 @@ public transportEquationProblem< Mesh, BoundaryCondition, RightHandSide, Differe
		typedef typename DifferentialOperator::VelocityFieldType VelocityFieldType;
		typedef SharedPointer< VelocityFieldType, DeviceType > VelocityFieldPointer;


		using typename BaseType::MeshType;
		using typename BaseType::MeshPointer;
		using typename BaseType::DofVectorType;
		@@ -47,6 +48,8 @@ public transportEquationProblem< Mesh, BoundaryCondition, RightHandSide, Differe
		using typename BaseType::MeshDependentDataType;
		using typename BaseType::MeshDependentDataPointer;

		//using BaseType::getExplicitRHS;

		static String getTypeStatic();

		String getPrologHeader() const;

examples/transport-equation/transportEquationProblemEoc_impl.h

+43 −138

Original line number	Diff line number	Diff line
		@@ -76,31 +76,43 @@ setup( const MeshPointer& meshPointer,
		typedef Functions::Analytic::Paraboloid< Dimensions, RealType > ParaboloidType;
		typedef Operators::Analytic::Heaviside< ParaboloidType > HeavisideParaboloidType;
		typedef Functions::OperatorFunction< HeavisideParaboloidType, ParaboloidType > InitialConditionType;
		String velocityFieldType = parameters.getParameter< String >( "velocity-field" );
		if( velocityFieldType == "constant" )
		{
		typedef Operators::Analytic::Shift< InitialConditionType > ShiftOperatorType;
		typedef Functions::OperatorFunction< ShiftOperatorType, InitialConditionType > ExactSolutionType;
		SharedPointer< ExactSolutionType, Devices::Host > exactSolution;
		if( ! exactSolution->setup( parameters, prefix ) )
		if( ! exactSolution->getFunction().setup( parameters, prefix ) )
		return false;
		/*Functions::MeshFunctionEvaluator< MeshFunction, ExactSolutionType > evaluator;
		Containers::StaticVector< Dimensions, RealType > velocity;
		for( int i = 0; i < Dimensions; i++ )
		velocity[ i ] = parameters.getParameter< double >( "velocity-field-" + String( i ) + "-constant" );

		Functions::MeshFunctionEvaluator< MeshFunction, ExactSolutionType > evaluator;
		RealType time( 0.0 );
		int step( 0 );
		exactSolution->getOperator().setShift( 0.0 * velocity );
		evaluator.evaluate( u, exactSolution, time );
		FileName fileName;
		fileName.setFileNameBase( "exact-u-" );
		fileName.setExtension( "tnl" );
		fileName.setIndex( step );
		if( ! u.save( fileName.getFileName() ) )
		if( ! u->save( fileName.getFileName() ) )
		return false;
		while( time < finalTime )
		{
		time += snapshotPeriod;
		if( time > finalTime )
		time = finalTime;
		exactSolution->getOperator().setShift( time * velocity );
		std::cerr << time * velocity << std::endl;
		std::cerr << exactSolution->getOperator().getShift() << std::endl;
		evaluator.evaluate( u, exactSolution, time );
		fileName.setIndex( ++step );
		if( ! u.save( fileName.getFileName() ) )
		if( ! u->save( fileName.getFileName() ) )
		return false;
		}*/
		}
		}
		}

		return true;
		@@ -132,124 +144,17 @@ setInitialCondition( const Config::ParameterContainer& parameters,
		MeshDependentDataPointer& meshDependentData )
		{
		this->bindDofs( meshPointer, dofs );
		const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" );
		if( ! this->uPointer->boundLoad( initialConditionFile ) )
		{
		std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << std::endl;
		return false;
		}
		return true;
		}

		#ifdef UNDEF
		template< typename Mesh,
		typename BoundaryCondition,
		typename RightHandSide,
		typename DifferentialOperator >
		template< typename Matrix >
		bool
		transportEquationProblemEoc< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
		setupLinearSystem( const MeshPointer& mesh,
		Matrix& matrix )
		{
		/*const IndexType dofs = this->getDofs( mesh );
		typedef typename Matrix::ObjectType::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
		SharedPointer< CompressedRowsLengthsVectorType > rowLengths;
		if( ! rowLengths->setSize( dofs ) )
		return false;
		Matrices::MatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, CompressedRowsLengthsVectorType > matrixSetter;
		matrixSetter.template getCompressedRowsLengths< typename Mesh::Cell >( mesh,
		differentialOperatorPointer,
		boundaryConditionPointer,
		rowLengths );
		matrix->setDimensions( dofs, dofs );
		if( ! matrix->setCompressedRowsLengths( *rowLengths ) )
		return false;*/
		return true;
		}

		template< typename Mesh,
		typename BoundaryCondition,
		typename RightHandSide,
		typename DifferentialOperator >
		bool
		transportEquationProblemEoc< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
		makeSnapshot( const RealType& time,
		const IndexType& step,
		const MeshPointer& mesh,
		DofVectorPointer& dofs,
		MeshDependentDataPointer& meshDependentData )
		{
		std::cout << std::endl << "Writing output at time " << time << " step " << step << "." << std::endl;
		this->bindDofs( mesh, dofs );
		//const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" );
		FileName fileName;
		fileName.setFileNameBase( "u-" );
		fileName.setFileNameBase( "exact-u-" );
		fileName.setExtension( "tnl" );
		fileName.setIndex( step );
		if( ! dofs->save( fileName.getFileName() ) )
		return false;
		return true;
		}

		template< typename Mesh,
		typename BoundaryCondition,
		typename RightHandSide,
		typename DifferentialOperator >
		void
		transportEquationProblemEoc< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
		getExplicitRHS( const RealType& time,
		const RealType& tau,
		const MeshPointer& mesh,
		DofVectorPointer& _u,
		DofVectorPointer& _fu,
		MeshDependentDataPointer& meshDependentData )
		fileName.setIndex( 0 );
		if( ! this->uPointer->boundLoad( fileName.getFileName() ) )
		{
		/****
		* If you use an explicit solver like Euler or Merson, you
		* need to implement this method. Compute the right-hand side of
		*
		* d/dt u(x) = fu( x, u )
		*
		* You may use supporting mesh dependent data if you need.
		*/
		typedef typename MeshType::Cell Cell;
		int count = ::sqrt(mesh->template getEntitiesCount< Cell >());
		this->bindDofs( mesh, _u );
		Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
		SharedPointer< MeshFunctionType > u( mesh, _u );
		SharedPointer< MeshFunctionType > fu( mesh, _fu );
		differentialOperatorPointer->setTau(tau);
		differentialOperatorPointer->setVelocityField( this->velocityField );
		explicitUpdater.template update< typename Mesh::Cell >( time,
		mesh,
		this->differentialOperatorPointer,
		this->boundaryConditionPointer,
		this->rightHandSidePointer,
		u,
		fu );
		/*BoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
		boundaryConditionsSetter.template apply< typename Mesh::Cell >(
		this->boundaryCondition,
		time + tau,
		u ); */
		std::cerr << "I am not able to load the initial condition from the file " << fileName.getFileName() << "." << std::endl;
		return false;
		}
		template< typename Mesh,
		typename BoundaryCondition,
		typename RightHandSide,
		typename DifferentialOperator >
		template< typename Matrix >
		void
		transportEquationProblemEoc< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
		assemblyLinearSystem( const RealType& time,
		const RealType& tau,
		const MeshPointer& mesh,
		DofVectorPointer& _u,
		Matrix& matrix,
		DofVectorPointer& b,
		MeshDependentDataPointer& meshDependentData )
		{
		return true;
		}

		#endif

		} // namespace TNL

src/TNL/Config/ParameterContainer.h

+1 −1

Original line number	Diff line number	Diff line
		@@ -63,7 +63,7 @@ class ParameterContainer

		template< class T > bool getParameter( const String& name,
		T& value,
		bool verbose = false ) const
		bool verbose = true ) const
		{
		int i;
		const int size = parameters. getSize();

src/TNL/Containers/StaticVector.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -23,6 +23,8 @@ class StaticVector : public Containers::StaticArray< Size, Real >
		typedef StaticVector< Size, Real > ThisType;
		enum { size = Size };

		using Containers::StaticArray< Size, Real >::operator=;

		__cuda_callable__
		StaticVector();