Commit 28b18c26 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Refactoring explicit updater and matrix assembler.

parent 5ddaa9bc
Loading
Loading
Loading
Loading
+12 −20
Original line number Diff line number Diff line
@@ -280,25 +280,17 @@ getExplicitUpdate( const RealType& time,
	       *(_u[i*count+j-1]-_u[(i-2)*count+j-1]))
            ;}
  }
   else if (this->velocityType == "advection")
*/  { 
   else if (this->velocityType == "advection")*/
   { 
      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); 
   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 ); */
      explicitUpdater.setDifferentialOperator( this->differentialOperatorPointer );
      explicitUpdater.setBoundaryConditions( this->boundaryConditionPointer );
      explicitUpdater.setRightHandSide( this->rightHandSidePointer );
      explicitUpdater.template update< typename Mesh::Cell >( time, tau, mesh, u, fu );
   }
}
template< typename Mesh,
+9 −14
Original line number Diff line number Diff line
@@ -297,13 +297,10 @@ getExplicitUpdate( const RealType& time,
   lF1DMomentum->setVelocity( *velocity);
   lF1DMomentum->setPressure( *pressure);
   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Momentum, BoundaryCondition, RightHandSide > explicitUpdaterMomentum;
   explicitUpdaterMomentum.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           lF1DMomentum,
                                                           this->boundaryConditionsPointer,
                                                           this->rightHandSidePointer,
                                                           uRhoVelocity,
                                                           fuRhoVelocity );
   explicitUpdaterMomentum.setDifferentialOperator( lF1DMomentum );
   explicitUpdaterMomentum.setBoundaryConditions( this->boundaryConditionsPointer );
   explicitUpdaterMomentum.setRightHandSide( this->rightHandSidePointer );
   explicitUpdaterMomentum.template update< typename Mesh::Cell >( time, tau, mesh, uRhoVelocity, fuRhoVelocity );
   
   std::cout << "explicitRHSenergy" << std::endl;
   //energy
@@ -311,13 +308,11 @@ getExplicitUpdate( const RealType& time,
   lF1DEnergy->setPressure( *pressure);
   lF1DEnergy->setVelocity( *velocity);
   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Energy, BoundaryCondition, RightHandSide > explicitUpdaterEnergy;
   explicitUpdaterEnergy.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           lF1DEnergy,
                                                           this->boundaryConditionsPointer,
                                                           this->rightHandSidePointer,
                                                           uEnergy,
                                                           fuEnergy );  
   explicitUpdaterEnergy.setDifferentialOperator( lF1DEnergy );
   explicitUpdaterEnergy.setBoundaryConditions( this->boundaryConditionsPointer );
   explicitUpdaterEnergy.setRightHandSide( this->rightHandSidePointer );

   explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, tau, mesh, uEnergy, fuEnergy );  
 
 }

+18 −36
Original line number Diff line number Diff line
@@ -280,13 +280,10 @@ getExplicitUpdate( const RealType& time,
   lF2DContinuity->setVelocityX( *velocityX );
   lF2DContinuity->setVelocityY( *velocityY );
   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Continuity, BoundaryCondition, RightHandSide > explicitUpdaterContinuity; 
   explicitUpdaterContinuity.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           lF2DContinuity,
                                                           this->boundaryConditionPointer,
                                                           this->rightHandSidePointer,
                                                           uRho,
                                                           fuRho );
   explicitUpdaterContinuity.setDifferentialOperator( lF2DContinuity );
   explicitUpdaterContinuity.setBoundaryConditions( this->boundaryConditionPointer );
   explicitUpdaterContinuity.setRightHandSide( this->rightHandSidePointer );
   explicitUpdaterContinuity.template update< typename Mesh::Cell >( time, tau, mesh, uRho, fuRho );

   //rhoVelocityX
   lF2DMomentumX->setTau(tau);
@@ -294,13 +291,10 @@ getExplicitUpdate( const RealType& time,
   lF2DMomentumX->setVelocityY( *velocityY );
   lF2DMomentumX->setPressure( *pressure );
   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumX, BoundaryCondition, RightHandSide > explicitUpdaterMomentumX; 
   explicitUpdaterMomentumX.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           lF2DMomentumX,
                                                           this->boundaryConditionPointer,
                                                           this->rightHandSidePointer,
                                                           uRhoVelocityX,
                                                           fuRhoVelocityX );
   explicitUpdaterMomentumX.setDifferentialOperator( lF2DMomentumX );
   explicitUpdaterMomentumX.setBoundaryConditions( this->boundaryConditionPointer );
   explicitUpdaterMomentumX.setRightHandSide( this->rightHandSidePointer );   
   explicitUpdaterMomentumX.template update< typename Mesh::Cell >( time, tau, mesh, uRhoVelocityX, fuRhoVelocityX );

   //rhoVelocityY
   lF2DMomentumY->setTau(tau);
@@ -308,13 +302,10 @@ getExplicitUpdate( const RealType& time,
   lF2DMomentumY->setVelocityY( *velocityY );
   lF2DMomentumY->setPressure( *pressure );
   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumY, BoundaryCondition, RightHandSide > explicitUpdaterMomentumY;
   explicitUpdaterMomentumY.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           lF2DMomentumY,
                                                           this->boundaryConditionPointer,
                                                           this->rightHandSidePointer,
                                                           uRhoVelocityY,
                                                           fuRhoVelocityY );
   explicitUpdaterMomentumY.setDifferentialOperator( lF2DMomentumY );
   explicitUpdaterMomentumY.setBoundaryConditions( this->boundaryConditionPointer );
   explicitUpdaterMomentumY.setRightHandSide( this->rightHandSidePointer );
   explicitUpdaterMomentumY.template update< typename Mesh::Cell >( time, tau, mesh, uRhoVelocityY, fuRhoVelocityY );
  
   //energy
   lF2DEnergy->setTau(tau);
@@ -322,20 +313,11 @@ getExplicitUpdate( const RealType& time,
   lF2DEnergy->setVelocityY( *velocityY ); 
   lF2DEnergy->setPressure( *pressure );
   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Energy, BoundaryCondition, RightHandSide > explicitUpdaterEnergy;
   explicitUpdaterEnergy.template update< typename Mesh::Cell >( time,
                                                           mesh,
                                                           lF2DEnergy,
                                                           this->boundaryConditionPointer,
                                                           this->rightHandSidePointer,
                                                           uEnergy,
                                                           fuEnergy );
   explicitUpdaterEnergy.setDifferentialOperator( lF2DEnergy );
   explicitUpdaterEnergy.setBoundaryConditions( this->boundaryConditionPointer );
   explicitUpdaterEnergy.setRightHandSide( this->rightHandSidePointer );
   explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, tau, mesh, uEnergy, fuEnergy );

/*
   BoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 
   boundaryConditionsSetter.template apply< typename Mesh::Cell >( 
      this->boundaryCondition, 
      time + tau, 
       u );*/
}

template< typename Mesh,
@@ -433,7 +415,7 @@ postIterate( const RealType& time,
   euler2DPressure.setRho(uRho);
//   OperatorFunction< euler2DPressure, MeshFunction, void, time > OFPressure;
//   pressure = OFPressure;

   return true;
}

} // namespace TNL
+10 −0
Original line number Diff line number Diff line
@@ -22,6 +22,8 @@
#include <TNL/Functions/MeshFunction.h>
#include <TNL/Timer.h>
#include <TNL/Solvers/PDE/ExplicitUpdater.h>
#include <TNL/Solvers/PDE/LinearSystemAssembler.h>
#include <TNL/Solvers/PDE/BackwardTimeDiscretisation.h>

namespace TNL {
namespace Problems {
@@ -120,6 +122,14 @@ class HeatEquationProblem : public PDEProblem< Mesh,
         Timer gpuTransferTimer;
         
         Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
         
         Solvers::PDE::LinearSystemAssembler< Mesh, 
                                              MeshFunctionType,
                                              DifferentialOperator,
                                              BoundaryCondition,
                                              RightHandSide,
                                              Solvers::PDE::BackwardTimeDiscretisation,
                                              DofVectorType > systemAssembler;
};

} // namespace Problems
+9 −60
Original line number Diff line number Diff line
@@ -21,8 +21,6 @@
#include <TNL/Matrices/MultidiagonalMatrixSetter.h>
#include <TNL/Logger.h>
#include <TNL/Solvers/PDE/BoundaryConditionsSetter.h>
#include <TNL/Solvers/PDE/LinearSystemAssembler.h>
#include <TNL/Solvers/PDE/BackwardTimeDiscretisation.h>

#include "HeatEquationProblem.h"

@@ -220,26 +218,10 @@ getExplicitUpdate( const RealType& time,
   
   this->bindDofs( meshPointer, uDofs );
   MeshFunctionPointer fuPointer( meshPointer, fuDofs );
   explicitUpdater.template update< typename Mesh::Cell >(
      time,
      meshPointer,
      this->differentialOperatorPointer,
      this->boundaryConditionPointer,
      this->rightHandSidePointer,
      this->uPointer,
      fuPointer );
   //std::cerr << "******************************************************************************************" << std::endl;
   //std::cerr << "******************************************************************************************" << std::endl;
   //std::cerr << "******************************************************************************************" << std::endl;
   /*Solvers::PDE::BoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
   boundaryConditionsSetter.template apply< typename Mesh::Cell >(
      this->boundaryConditionPointer,
      time + tau,
      this->uPointer );*/
   
   /*uPointer->write( "u.txt", "gnuplot" );
   fuPointer->write( "fu.txt", "gnuplot" );
   getchar();*/
   this->explicitUpdater.setDifferentialOperator( this->differentialOperatorPointer ),
   this->explicitUpdater.setBoundaryConditions( this->boundaryConditionPointer ),
   this->explicitUpdater.setRightHandSide( this->rightHandSidePointer ),
   this->explicitUpdater.template update< typename Mesh::Cell >( time, tau, meshPointer, this->uPointer, fuPointer );
}

template< typename Mesh,
@@ -258,49 +240,16 @@ assemblyLinearSystem( const RealType& time,
                      MeshDependentDataPointer& meshDependentData )
{
   this->bindDofs( meshPointer, dofsPointer );
   Solvers::PDE::LinearSystemAssembler< Mesh,
                             MeshFunctionType,
                             DifferentialOperator,
                             BoundaryCondition,
                             RightHandSide,
                             Solvers::PDE::BackwardTimeDiscretisation,
                             typename MatrixPointer::ObjectType,
                             DofVectorType > systemAssembler;
   systemAssembler.template assembly< typename Mesh::Cell >(
   this->systemAssembler.setDifferentialOperator( this->differentialOperatorPointer );
   this->systemAssembler.setBoundaryConditions( this->boundaryConditionPointer );
   this->systemAssembler.setRightHandSide( this->rightHandSidePointer );
   this->systemAssembler.template assembly< typename Mesh::Cell, typename MatrixPointer::ObjectType >( 
      time,
      tau,
      meshPointer,
      this->differentialOperatorPointer,
      this->boundaryConditionPointer,
      this->rightHandSidePointer,
      this->uPointer,
      matrixPointer,
      bPointer );
   //matrixPointer->print( std::cout );
   //getchar();
   /*cout << endl << b << endl;
   cout << endl << u << endl;
   abort();*/
   /*cout << "Matrix multiplication test ..." << std::endl;
   Vector< RealType, DeviceType, IndexType > y;
   y.setLike( u );
   Timer timer;
   timer.reset();
   timer.start();
   for( int i = 0; i < 100; i++ )
      matrix.vectorProduct( u, y );
   timer.stop();
  std::cout << "The time is " << timer.getRealTime();
  std::cout << "Scalar product test ..." << std::endl;
   timer.reset();
   RealType a;
   timer.start();
   for( int i = 0; i < 100; i++ )
      a = y.scalarProduct( u );
   timer.stop();
  std::cout << "The time is " << timer.getRealTime();
  std::cout << std::endl;
   abort();*/
}

} // namespace Problems
Loading