Implementing exact gradient norm. (f5d084b5) · Commits · TNL / tnl-dev

examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h

+2 −2

Original line number	Diff line number	Diff line
		@@ -33,7 +33,7 @@
		#include <operators/diffusion/tnlExactNonlinearDiffusion.h>
		#include <operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h>
		#include <operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h>
		#include <operators/operator-Q/tnlExactOperatorQ.h>
		#include <operators/geometric/tnlExactGradientNorm.h>

		//typedef tnlDefaultConfigTag BuildConfig;
		typedef tnlFastBuildConfig BuildConfig;
		@@ -78,7 +78,7 @@ class meanCurvatureFlowEocSetter
		typedef tnlFiniteVolumeOperatorQ<MeshType, Real, Index, 0> OperatorQ;
		typedef tnlFiniteVolumeNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;
		typedef tnlNonlinearDiffusion< MeshType, NonlinearOperator, Real, Index > ApproximateOperator;
		typedef tnlExactNonlinearDiffusion< tnlExactOperatorQ<Dimensions>, Dimensions > ExactOperator;
		typedef tnlExactNonlinearDiffusion< tnlExactGradientNorm< Dimensions >, Dimensions > ExactOperator;
		typedef tnlTestFunction< MeshType::meshDimensions, Real, Device > TestFunction;
		typedef tnlMeanCurvatureFlowEocRhs< ExactOperator, TestFunction, Dimensions > RightHandSide;
		typedef tnlStaticVector < MeshType::meshDimensions, Real > Vertex;

src/operators/fdm/tnlForwardFiniteDifference.h

+3 −0

Original line number	Diff line number	Diff line
		@@ -18,6 +18,9 @@
		#ifndef TNLFORWARDFINITEDIFFERENCE_H
		#define TNLFORWARDFINITEDIFFERENCE_H

		#include<mesh/tnlGrid.h>
		#include<operators/fdm/tnlFiniteDifferences.h>

		template< typename Mesh,
		int Xdifference = 0,
		int YDifference = 0,

src/operators/geometric/tnlExactGradientNorm.h

+83 −18

Original line number	Diff line number	Diff line
		@@ -65,8 +65,7 @@ class tnlExactGradientNorm< 1, Real >
		template< typename Function,
		int XDerivative = 0,
		int YDerivative = 0,
		int ZDerivative = 0,
		typename Real = typename Function::RealType >
		int ZDerivative = 0 >
		__cuda_callable__
		typename Function::RealType
		getPartialDerivative( const Function& function,
		@@ -76,12 +75,8 @@ class tnlExactGradientNorm< 1, Real >
		static_assert( XDerivative >= 0 && YDerivative >= 0 && ZDerivative >= 0,
		"Partial derivative must be non-negative integer." );
		static_assert( XDerivative < 2, "Partial derivative of higher order then 1 are not implemented yet." );
		typedef typename Function::RealType ReealType;
		typedef typename Function::RealType RealType;

		if( YDerivative != 0 \|\| ZDerivative != 0 )
		return 0.0;
		if( XDerivative == 0 )
		return this->operator()( function, vertex, time );
		if( XDerivative == 1 )
		{
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		@@ -89,6 +84,10 @@ class tnlExactGradientNorm< 1, Real >
		const RealType Q = sqrt( this->epsilonSquare + f_x * f_x );
		return ( f_x * f_xx ) / Q;
		}
		if( XDerivative == 0 )
		return this->operator()( function, v, time );
		if( YDerivative != 0 \|\| ZDerivative != 0 )
		return 0.0;
		}

		protected:
		@@ -126,22 +125,47 @@ class tnlExactGradientNorm< 2, Real >
		const typename Function::VertexType& v,
		const typename Function::RealType& time = 0.0 ) const
		{

		typedef typename Function::RealType RealType;
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		return sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y );
		}

		template< typename Function,
		int XDerivative = 0,
		int YDerivative = 0,
		int ZDerivative = 0,
		typename Real = typename Function::RealType >
		int ZDerivative = 0 >
		__cuda_callable__
		typename Function::RealType
		getPartialDerivative( const Function& function,
		const typename Function::Vertex& v,
		const Real& time = 0.0,
		const Real& eps = 1.0 ) const
		const Real& time = 0.0 ) const
		{
		static_assert( XDerivative >= 0 && YDerivative >= 0 && ZDerivative >= 0,
		"Partial derivative must be non-negative integer." );
		static_assert( XDerivative < 2 && YDerivative < 2, "Partial derivative of higher order then 1 are not implemented yet." );
		typedef typename Function::RealType RealType;

		if( XDerivative == 1 && YDerivative == 0 )
		{
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		const RealType f_xx = function.template getPartialDerivative< 2, 0, 0 >( v, time );
		const RealType f_xy = function.template getPartialDerivative< 1, 1, 0 >( v, time );
		return ( f_x * f_xx + f_y * f_xy ) / sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y );
		}
		if( XDerivative == 0 && YDerivative == 1 )
		{
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		const RealType f_xy = function.template getPartialDerivative< 1, 1, 0 >( v, time );
		const RealType f_yy = function.template getPartialDerivative< 0, 2, 0 >( v, time );
		return ( f_x * f_xy + f_y * f_yy ) / sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y );
		}
		if( XDerivative == 0 && YDerivative == 0 )
		return this->operator()( function, v, time );
		if( ZDerivative > 0 )
		return 0.0;
		}

		protected:
		@@ -175,19 +199,64 @@ class tnlExactGradientNorm< 3, Real >
		const typename Function::VertexType& v,
		const typename Function::RealType& time = 0.0 ) const
		{
		typedef typename Function::RealType RealType;
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		const RealType f_z = function.template getPartialDerivative< 0, 0, 1 >( v, time );
		return sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y + f_z * f_z );
		}

		template< typename Function,
		int XDerivative = 0,
		int YDerivative = 0,
		int ZDerivative = 0,
		typename Real = typename Function::RealType >
		int ZDerivative = 0 >
		__cuda_callable__
		typename Function::RealType
		getPartialDerivative( const Function& function,
		const typename Function::Vertex& v,
		const Real& time = 0.0 ) const
		{
		static_assert( XDerivative >= 0 && YDerivative >= 0 && ZDerivative >= 0,
		"Partial derivative must be non-negative integer." );
		static_assert( XDerivative < 2 && YDerivative < 2 && ZDerivative < 2,
		"Partial derivative of higher order then 1 are not implemented yet." );
		typedef typename Function::RealType RealType;

		if( XDerivative == 1 && YDerivative == 0 && ZDerivative == 0 )
		{
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		const RealType f_z = function.template getPartialDerivative< 0, 0, 1 >( v, time );
		const RealType f_xx = function.template getPartialDerivative< 2, 0, 0 >( v, time );
		const RealType f_xy = function.template getPartialDerivative< 1, 1, 0 >( v, time );
		const RealType f_xz = function.template getPartialDerivative< 1, 0, 1 >( v, time );
		return ( f_x * f_xx + f_y * f_xy + f_z * f_xz ) /
		sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y + f_z * f_z );
		}
		if( XDerivative == 0 && YDerivative == 1 && ZDerivative == 0 )
		{
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		const RealType f_z = function.template getPartialDerivative< 0, 0, 1 >( v, time );
		const RealType f_xy = function.template getPartialDerivative< 1, 1, 0 >( v, time );
		const RealType f_yy = function.template getPartialDerivative< 0, 2, 0 >( v, time );
		const RealType f_yz = function.template getPartialDerivative< 0, 1, 1 >( v, time );
		return ( f_x * f_xy + f_y * f_yy + f_z * f_yz ) /
		sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y + f_z * f_z );
		}
		if( XDerivative == 0 && YDerivative == 0 && ZDerivative == 1 )
		{
		const RealType f_x = function.template getPartialDerivative< 1, 0, 0 >( v, time );
		const RealType f_y = function.template getPartialDerivative< 0, 1, 0 >( v, time );
		const RealType f_z = function.template getPartialDerivative< 0, 0, 1 >( v, time );
		const RealType f_xz = function.template getPartialDerivative< 1, 0, 1 >( v, time );
		const RealType f_yz = function.template getPartialDerivative< 0, 1, 1 >( v, time );
		const RealType f_zz = function.template getPartialDerivative< 0, 0, 2 >( v, time );
		return ( f_x * f_xz + f_y * f_yz + f_z * f_zz ) /
		sqrt( this->epsilonSquare + f_x * f_x + f_y * f_y + f_z * f_z );
		}
		if( XDerivative == 0 && YDerivative == 0 && ZDerivative == 0 )
		return this->operator()( function, v, time );
		}

		protected:
		@@ -195,8 +264,4 @@ class tnlExactGradientNorm< 3, Real >
		Real epsilonSquare;
		};


		#include <operators/operator-Q/tnlExactGradientNorm_impl.h>


		#endif /* TNLEXACTGRADIENTTNORM_H */

src/operators/geometric/tnlFDMGradientNorm.h

+8 −3

Original line number	Diff line number	Diff line
		@@ -18,6 +18,8 @@
		#ifndef TNLFDMGRADIENTNORM_H
		#define TNLFDMGRADIENTNORM_H

		#include<operators/fdm/tnlForwardFiniteDifference.h>

		template< typename Mesh,
		template< typename, int, int, int, typename, typename > class DifferenceOperatorTemplate = tnlForwardFiniteDifference,
		typename Real = typename Mesh::RealType,
		@@ -32,7 +34,8 @@ template< typename MeshReal,
		template< typename, int, int, int, typename, typename > class DifferenceOperatorTemplate,
		typename Real,
		typename Index >
		class tnlFDMGradientNorm< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
		class tnlFDMGradientNorm< tnlGrid< 1,MeshReal, Device, MeshIndex >, DifferenceOperatorTemplate, Real, Index >
		: public tnlDomain< 1, MeshInteriorDomain >
		{
		public:

		@@ -81,7 +84,8 @@ template< typename MeshReal,
		template< typename, int, int, int, typename, typename > class DifferenceOperatorTemplate,
		typename Real,
		typename Index >
		class tnlFDMGradientNorm< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
		class tnlFDMGradientNorm< tnlGrid< 2,MeshReal, Device, MeshIndex >, DifferenceOperatorTemplate, Real, Index >
		: public tnlDomain< 2, MeshInteriorDomain >
		{
		public:

		@@ -136,7 +140,8 @@ template< typename MeshReal,
		template< typename, int, int, int, typename, typename > class DifferenceOperatorTemplate,
		typename Real,
		typename Index >
		class tnlFDMGradientNorm< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
		class tnlFDMGradientNorm< tnlGrid< 3, MeshReal, Device, MeshIndex >, DifferenceOperatorTemplate, Real, Index >
		: public tnlDomain< 3, MeshInteriorDomain >
		{
		public:

tests/unit-tests/operators/geometric/tnlFDMGradientNormTest.h

+94 −12

Original line number	Diff line number	Diff line
		@@ -20,13 +20,92 @@

		#include <operators/geometric/tnlFDMGradientNorm.h>
		#include <operators/geometric/tnlExactGradientNorm.h>
		#include <operators/fdm/tnlBackwardFiniteDifference.h>
		#include <operators/fdm/tnlForwardFiniteDifference.h>
		#include <operators/fdm/tnlCentralFiniteDifference.h>
		#include "../../tnlUnitTestStarter.h"
		#include "../tnlPDEOperatorEocTester.h"

		template< int Dimensions,
		typename Real,
		typename Device,
		typename Index,
		typename TestFunction >
		class tnlPDEOperatorEocTestResult<
		tnlFDMGradientNorm< tnlGrid< Dimensions, Real, Device, Index >,
		tnlForwardFiniteDifference,
		Real,
		Index >,
		TestFunction >
		{
		public:
		static Real getL1Eoc() { return ( Real ) 1.0; };
		static Real getL1Tolerance() { return ( Real ) 0.05; };

		static Real getL2Eoc() { return ( Real ) 1.0; };
		static Real getL2Tolerance() { return ( Real ) 0.05; };

		static Real getMaxEoc() { return ( Real ) 1.0; };
		static Real getMaxTolerance() { return ( Real ) 0.05; };

		};

		template< int Dimensions,
		typename Real,
		typename Device,
		typename Index,
		typename TestFunction >
		class tnlPDEOperatorEocTestResult<
		tnlFDMGradientNorm< tnlGrid< Dimensions, Real, Device, Index >,
		tnlBackwardFiniteDifference,
		Real,
		Index >,
		TestFunction >
		{
		public:
		static Real getL1Eoc() { return ( Real ) 1.0; };
		static Real getL1Tolerance() { return ( Real ) 0.05; };

		static Real getL2Eoc() { return ( Real ) 1.0; };
		static Real getL2Tolerance() { return ( Real ) 0.05; };

		static Real getMaxEoc() { return ( Real ) 1.0; };
		static Real getMaxTolerance() { return ( Real ) 0.05; };

		};

		template< int Dimensions,
		typename Real,
		typename Device,
		typename Index,
		typename TestFunction >
		class tnlPDEOperatorEocTestResult<
		tnlFDMGradientNorm< tnlGrid< Dimensions, Real, Device, Index >,
		tnlCentralFiniteDifference,
		Real,
		Index >,
		TestFunction >
		{
		public:
		static Real getL1Eoc() { return ( Real ) 2.0; };
		static Real getL1Tolerance() { return ( Real ) 0.05; };

		static Real getL2Eoc() { return ( Real ) 2.0; };
		static Real getL2Tolerance() { return ( Real ) 0.05; };

		static Real getMaxEoc() { return ( Real ) 2.0; };
		static Real getMaxTolerance() { return ( Real ) 0.05; };

		};



		template< typename Mesh,
		typename Function,
		typename Operator,
		int MeshSize,
		bool Verbose >
		void testDifferenceOperator()
		bool testDifferenceOperator()
		{
		typedef tnlExactGradientNorm< Mesh::meshDimensions > ExactOperator;
		return tnlUnitTestStarter::run<
		@@ -43,30 +122,33 @@ template< typename Mesh,
		typename Function,
		int MeshSize,
		bool Verbose >
		void setDifferenceOperator()
		bool setDifferenceOperator()
		{
		typedef tnlFDMGradientNorm< Mesh, tnlForwardFiniteDifference > ForwardGradientNorm;
		typedef tnlFDMGradientNorm< Mesh, tnlBackwardFiniteDifference > BackwardGradientNorm;
		typedef tnlFDMGradientNorm< Mesh, tnlCentralFiniteDifference > CentralGradientNorm;
		return ( testDifferenceOperator< Mesh, Function, ForwardGradientNorm, MeshSize, Verboze >() &&
		testDifferenceOperator< Mesh, Function, BackwardGradientNorm, MeshSize, Verboze >() &&
		testDifferenceOperator< Mesh, Function, CentralGradientNorm, MeshSize, Verboze >() );
		return ( testDifferenceOperator< Mesh, Function, ForwardGradientNorm, MeshSize, Verbose >() &&
		testDifferenceOperator< Mesh, Function, BackwardGradientNorm, MeshSize, Verbose >() &&
		testDifferenceOperator< Mesh, Function, CentralGradientNorm, MeshSize, Verbose >() );
		}

		template< typename Mesh,
		int MeshSize,
		bool Verbose >
		void setFunction()
		bool setFunction()
		{
		const int Dimensions = Mesh::meshDimensions;
		typedef tnlExpBumpFunction< Dimensions, RealType > Function;
		return setDiferenceOperator< Mesh, Function, MeshSize, Verbose >();
		typedef tnlExpBumpFunction< Dimensions, double > Function;
		return setDifferenceOperator< Mesh, Function, MeshSize, Verbose >();
		}

		template< int MeshSize,
		template< typename Device,
		int MeshSize,
		bool Verbose >
		void setGrid()
		bool setGrid()
		{
		typedef double MeshReal;
		typedef int MeshIndex;
		typedef tnlGrid< 1, MeshReal, Device, MeshIndex > Grid1D;
		typedef tnlGrid< 2, MeshReal, Device, MeshIndex > Grid2D;
		typedef tnlGrid< 3, MeshReal, Device, MeshIndex > Grid3D;
		@@ -80,7 +162,7 @@ int main( int argc, char* argv[] )
		const int meshSize( 32 );
		const bool verbose( false );
		#ifdef HAVE_CPPUNIT
		return setGrid< MeshSize, Verbose >();
		return setGrid< tnlHost, meshSize, verbose >();
		#else
		return EXIT_FAILURE;
		#endif