residualbased_eulerian_convdiff_strategy.h

Go to the documentation of this file.

// KRATOS ___ ___  _  ___   __   ___ ___ ___ ___
//       / __/ _ \| \| \ \ / /__|   \_ _| __| __|
//      | (_| (_) | .` |\ V /___| |) | || _|| _|
//       \___\___/|_|\_| \_/    |___/___|_| |_|  APPLICATION
//
//  License: BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:  Riccardo Rossi
//
 
#if !defined(KRATOS_RESIDUALBASED_EULERIAN_CONVECTION_DIFFUSION_STRATEGY )
#define  KRATOS_RESIDUALBASED_EULERIAN_CONVECTION_DIFFUSION_STRATEGY
 
 
/* System includes */
 
 
/* External includes */
 
/* Project includes */
#include "includes/define.h"
#include "containers/model.h"
#include "includes/model_part.h"
#include "solving_strategies/strategies/implicit_solving_strategy.h"
#include "solving_strategies/strategies/residualbased_linear_strategy.h"
#include "solving_strategies/schemes/residualbased_incrementalupdate_static_scheme.h"
#include "convection_diffusion_application.h"
#include "solving_strategies/builder_and_solvers/residualbased_elimination_builder_and_solver_componentwise.h"
#include "includes/convection_diffusion_settings.h"
//#include "custom_utilities/convection_diffusion_settings.h"
 
 
 
namespace Kratos
{
 
 
template<class TSparseSpace,
         class TDenseSpace,
         class TLinearSolver
         >
class ResidualBasedEulerianConvectionDiffusionStrategy
    : public ImplicitSolvingStrategy<TSparseSpace,TDenseSpace,TLinearSolver>
{
public:
    KRATOS_CLASS_POINTER_DEFINITION( ResidualBasedEulerianConvectionDiffusionStrategy );
 
    typedef ImplicitSolvingStrategy<TSparseSpace,TDenseSpace,TLinearSolver> BaseType;
 
    typedef typename BaseType::TDataType TDataType;
 
    //typedef typename BaseType::DofSetType DofSetType;
 
    typedef typename BaseType::DofsArrayType DofsArrayType;
 
    typedef typename BaseType::TSystemMatrixType TSystemMatrixType;
 
    typedef typename BaseType::TSystemVectorType TSystemVectorType;
 
    typedef typename BaseType::LocalSystemVectorType LocalSystemVectorType;
 
    typedef typename BaseType::LocalSystemMatrixType LocalSystemMatrixType;
 
 
 
    ResidualBasedEulerianConvectionDiffusionStrategy(
        ModelPart& model_part,
        typename TLinearSolver::Pointer pNewLinearSolver,
        bool ReformDofAtEachIteration = false,
        int dimension = 3
    )
        :
        ImplicitSolvingStrategy<TSparseSpace,TDenseSpace,TLinearSolver>(model_part,false),
        mrModelPart(model_part)
    {
        KRATOS_TRY
 
        GenerateMeshPart(dimension);
        mdimension = dimension;
        mOldDt = 0.00;
 
        //ProcessInfo& rCurrentProcessInfo = BaseType::GetModelPart().GetProcessInfo();
        //ConvectionDiffusionSettings::Pointer my_settings = rCurrentProcessInfo.GetValue(CONVECTION_DIFFUSION_SETTINGS);
        //(mpConvectionModelPart->GetProcessInfo()).GetValue(CONVECTION_DIFFUSION_SETTINGS) = my_settings;
        Check();
 
 
        //initializing fractional velocity solution step
        typedef Scheme< TSparseSpace,  TDenseSpace > SchemeType;
        typename SchemeType::Pointer pscheme = typename SchemeType::Pointer
                                               ( new ResidualBasedIncrementalUpdateStaticScheme< TSparseSpace,  TDenseSpace >() );
 
        bool CalculateReactions = false;
        bool CalculateNormDxFlag = true;
 
        //choosing the solving strategy
//          mstep1 = typename BaseType::Pointer(
//              new ResidualBasedLinearStrategy<TSparseSpace,  TDenseSpace, TLinearSolver >
//              (model_part,pscheme,pNewLinearSolver,CalculateReactions,ReformDofAtEachIteration,CalculateNormDxFlag)  );
//          mstep1->SetEchoLevel(2);
 
 
        typedef typename BuilderAndSolver<TSparseSpace,TDenseSpace,TLinearSolver>::Pointer BuilderSolverTypePointer;
 
        //const Variable<double>& rUnknownVar= my_settings->GetUnknownVariable();
        //BuilderSolverTypePointer componentwise_build = BuilderSolverTypePointer(new   ResidualBasedEliminationBuilderAndSolverComponentwise<TSparseSpace,TDenseSpace,TLinearSolver,Variable<double> > (pNewLinearSolver,rUnknownVar) );
        //mstep1 = typename BaseType::Pointer( new ResidualBasedLinearStrategy<TSparseSpace,  TDenseSpace, TLinearSolver >              (*mpConvectionModelPart,pscheme,pNewLinearSolver,componentwise_build,CalculateReactions,ReformDofAtEachIteration,CalculateNormDxFlag)  );
 
        BuilderSolverTypePointer pBuilderSolver = BuilderSolverTypePointer(new ResidualBasedBlockBuilderAndSolver<TSparseSpace,TDenseSpace,TLinearSolver>(pNewLinearSolver) );
        mstep1 = typename BaseType::Pointer( new ResidualBasedLinearStrategy<TSparseSpace,TDenseSpace,TLinearSolver >(
            *mpConvectionModelPart,
            pscheme,
            pBuilderSolver,
            CalculateReactions,
            ReformDofAtEachIteration,
            CalculateNormDxFlag));
 
        mstep1->SetEchoLevel(2);
 
        KRATOS_CATCH("")
    }
 
 
 
    virtual ~ResidualBasedEulerianConvectionDiffusionStrategy()
    {
        Model& current_model = mrModelPart.GetModel();
        current_model.DeleteModelPart("ConvDiffPart");
    }
 
    //*********************************************************************************
    //**********************************************************************
    double Solve() override
    {
      KRATOS_TRY
 
        //calculate the BDF coefficients
        //careful. Using information from the currentprocessinfo of the original model part (not the convection model part)
      //ProcessInfo& rCurrentProcessInfo = BaseType::GetModelPart().GetProcessInfo();
      //double Dt = rCurrentProcessInfo[DELTA_TIME];
      //int stationary= rCurrentProcessInfo[STATIONARY];
 
      //SOLVING THE PROBLEM
      double Dp_norm = mstep1->Solve();
 
      return Dp_norm;
      KRATOS_CATCH("")
    }
 
 
 
    void SetEchoLevel(int Level) override
    {
        mstep1->SetEchoLevel(Level);
    }
 
    void Clear() override
    {
        mstep1->Clear();
    }
 
    int Check() override
    {
        KRATOS_TRY
        ProcessInfo& rCurrentProcessInfo = BaseType::GetModelPart().GetProcessInfo();
        if (rCurrentProcessInfo.Has(CONVECTION_DIFFUSION_SETTINGS)==false)
            KRATOS_ERROR << "no CONVECTION_DIFFUSION_SETTINGS in model_part" << std::endl;
        //std::cout << "ConvDiff::Check(). If crashes, check CONVECTION_DIFFUSION_SETTINGS is defined" << std::endl;
 
        ConvectionDiffusionSettings::Pointer my_settings = rCurrentProcessInfo.GetValue(CONVECTION_DIFFUSION_SETTINGS);
 
        //DENSITY VARIABLE
        if(my_settings->IsDefinedDensityVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetDensityVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: Density Variable defined but not contained in the model part" << std::endl;
        }
        else
            std::cout << "No density variable assigned for ConvDiff. Assuming density=1" << std::endl;
 
        //DIFFUSION VARIABLE
        if(my_settings->IsDefinedDiffusionVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetDiffusionVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: Diffusion Variable defined but not contained in the model part" << std::endl;
        }
        else
            std::cout << "No diffusion variable assigned for ConvDiff. Assuming diffusivity=0" << std::endl;
 
        //UNKNOWN VARIABLE
        if(my_settings->IsDefinedUnknownVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetUnknownVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: Unknown Variable defined but not contained in the model part" << std::endl;
        }
        else
            KRATOS_ERROR << "ConvDiffSettings: Unknown Variable not defined!" << std::endl;
 
        //VOLUME SOURCE VARIABLE
        if(my_settings->IsDefinedVolumeSourceVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetVolumeSourceVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: VolumeSource Variable defined but not contained in the model part" << std::endl;
        }
        else
            std::cout << "No VolumeSource variable assigned for ConvDiff. Assuming VolumeSource=0" << std::endl;
 
 
        //SURFACE SOURCE VARIABLE
        //if(my_settings->IsDefinedSurfaceSourceVariable()==true)
        //{
        //  if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetSurfaceSourceVariable()) == false)
        //      KRATOS_ERROR << "ConvDiffSettings: SurfaceSource Variable defined but not contained in the model part" << std::endl;
        //}
        //else
        //  std::cout << "No SurfaceSource variable assigned for ConvDiff. Assuming SurfaceSource=0" << std::endl;
        if(my_settings->IsDefinedSurfaceSourceVariable()==true)
            KRATOS_ERROR << "ConvDiffSettings: SurfaceSource not yet implemented" << std::endl;
 
        //PROJECTION VARIABLE
        //if(my_settings->IsDefinedProjectionVariable()==true)
        //{
        //  if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetProjectionVariable()) == false)
        //      KRATOS_THROW_ERROR(std::logic_error, "ConvDiffSettings: Projection Variable defined but not contained in the model part", "");
        //}
        //else
        //  std::cout << "No Projection variable assigned for ConvDiff. Assuming Projection=0" << std::endl;
        if(my_settings->IsDefinedProjectionVariable()==true)
            KRATOS_ERROR << "ConvDiffSettings: ProjectionVariable not useed. Remove it" << std::endl;
 
        //CONVECTION VELOCITY VARIABLE
        //CURRENTLY WE ARE USING (VELOCITY -MESH_VELOCITY) TO CONVECT, so the ConvectionVariable must not be used:
        //if(my_settings->IsDefinedConvectionVariable()==true)
        //{
        //  if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetConvectionVariable()) == false)
        //      KRATOS_THROW_ERROR(std::logic_error, "ConvDiffSettings: Convection Variable defined but not contained in the model part", "");
        //}
        //else
        //  std::cout << "No Projection variable assigned for ConvDiff. Assuming Convection=0" << std::endl;
        if(my_settings->IsDefinedConvectionVariable()==true)
            KRATOS_ERROR << "ConvDiffSettings: ConvectionVariable not used. Use VelocityVariable instead" << std::endl;
 
        //MESH VELOCITY VARIABLE
        if(my_settings->IsDefinedMeshVelocityVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetMeshVelocityVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: MeshVelocity Variable defined but not contained in the model part" << std::endl;
        }
        else
            std::cout << "No MeshVelocity variable assigned for ConvDiff. Assuming MeshVelocity=0" << std::endl;
 
        //VELOCITY VARIABLE
        if(my_settings->IsDefinedVelocityVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetVelocityVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: Velocity Variable defined but not contained in the model part" << std::endl;
        }
        else
            std::cout << "No Velocity variable assigned for ConvDiff. Assuming Velocity=0" << std::endl;
 
        //TRANSFER COEFFICIENT VARIABLE
        //if(my_settings->IsDefinedTransferCoefficientVariable()==true)
        //{
        //  if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetTransferCoefficientVariable()) == false)
        //      KRATOS_THROW_ERROR(std::logic_error, "ConvDiffSettings: TransferCoefficient Variable defined but not contained in the model part", "");
        //}
        //else
        //  std::cout << "No TransferCoefficient variable assigned for ConvDiff. Assuming TransferCoefficient=0" << std::endl;
        if(my_settings->IsDefinedTransferCoefficientVariable()==true)
            KRATOS_ERROR << "ConvDiffSettings: TransferCoefficient not yet implemented" << std::endl;
 
        //SPECIFIC HEAT VARIABLE
        if(my_settings->IsDefinedSpecificHeatVariable()==true)
        {
            if (BaseType::GetModelPart().NodesBegin()->SolutionStepsDataHas(my_settings->GetSpecificHeatVariable()) == false)
                KRATOS_ERROR << "ConvDiffSettings: SpecificHeat Variable defined but not contained in the model part" << std::endl;
        }
        else
            std::cout << "No SpecificHeat variable assigned for ConvDiff. Assuming SpecificHeat=1" << std::endl;
 
        return 0;
 
        KRATOS_CATCH("")
 
    }
 
protected:
  virtual void GenerateMeshPart(int dimension)
  {
    Model& current_model = mrModelPart.GetModel();
    if(current_model.HasModelPart("ConvDiffPart"))
        current_model.DeleteModelPart("ConvDiffPart");
 
    // Generate
    mpConvectionModelPart = &(current_model.CreateModelPart("ConvDiffPart"));
 
 
    mpConvectionModelPart->SetProcessInfo(  BaseType::GetModelPart().pGetProcessInfo() );
    mpConvectionModelPart->SetBufferSize( BaseType::GetModelPart().GetBufferSize());
 
    //initializing mesh nodes
    mpConvectionModelPart->Nodes() = BaseType::GetModelPart().Nodes();
 
    //creating mesh elements
    ModelPart::ElementsContainerType& MeshElems = mpConvectionModelPart->Elements();
    Element::Pointer pElem;
 
    if(dimension == 2)
    {
        for(ModelPart::ElementsContainerType::iterator it= BaseType::GetModelPart().ElementsBegin(); it != BaseType::GetModelPart().ElementsEnd(); ++it)
         {
            Element::GeometryType& rGeom = it->GetGeometry();
            const unsigned int& NumNodes = rGeom.size();
 
            if (NumNodes == 3)
            {
              pElem = Element::Pointer(new EulerianConvectionDiffusionElement<2,3>(
                      (*it).Id(),
                      (*it).pGetGeometry(),
                      (*it).pGetProperties() ) );
              MeshElems.push_back(pElem);
            }
            else if(NumNodes == 4)
            {
                pElem = Element::Pointer(new EulerianConvectionDiffusionElement<2,4>(
                      (*it).Id(),
                      (*it).pGetGeometry(),
                      (*it).pGetProperties() ) );
              MeshElems.push_back(pElem);
            }
        }
    }
    else
    {
        for(ModelPart::ElementsContainerType::iterator it= BaseType::GetModelPart().ElementsBegin(); it != BaseType::GetModelPart().ElementsEnd(); ++it)
        {
            Element::GeometryType& rGeom = it->GetGeometry();
            const unsigned int& NumNodes = rGeom.size();
 
            if (NumNodes == 4)
            {
              pElem = Element::Pointer(new EulerianConvectionDiffusionElement<3,4>(
                      (*it).Id(),
                      (*it).pGetGeometry(),
                      (*it).pGetProperties() ) );
              MeshElems.push_back(pElem);
            }
            else if(NumNodes == 8)
            {
                pElem = Element::Pointer(new EulerianConvectionDiffusionElement<3,8>(
                      (*it).Id(),
                      (*it).pGetGeometry(),
                      (*it).pGetProperties() ) );
              MeshElems.push_back(pElem);
            }
         }
    }
  }
 
private:
    ModelPart& mrModelPart;
    ModelPart* mpConvectionModelPart;
    typename BaseType::Pointer mstep1;
    double mOldDt;
    int mdimension;
 
 
    ResidualBasedEulerianConvectionDiffusionStrategy(const ResidualBasedEulerianConvectionDiffusionStrategy& Other);
 
 
}; /* Class ResidualBasedEulerianConvectionDiffusionStrategy */
 
}  /* namespace Kratos.*/
 
#endif /* KRATOS_RESIDUALBASED_EULERIAN_CONVECTION_DIFFUSION_STRATEGY  defined */