How to Create Multi-step Elements

From KratosWiki
Jump to: navigation, search

On some occasions you want to solve several consecutive equations using the same mesh or model in Kratos, e.g. the following scheme:

u = f(x,y)  ---------->  v = g(u,x) ---------------> w = h(u,v)
  Step 1                   Step 2                      Step 3            

Suppose that in this scheme you want to solve a equation to get the value of u and then using this value to solve another equation and to obtain the value of v and finally compute the value of w using the values of u and v, therefore you need three steps to solve three equations and you need to know how to create this kind of multi-steps elements in Kratos. The first step is to modify the common custom element files in the folder of custom_elements. The first modification is to add three functions to introduce three steps for this element. We name our element as Element2D which you can change it according to the name of your application. The first modification is in the function of CalculateLocalSystem in Element2D.cpp:

Contents

C++ files

Defining LocalSystem function

Element cpp file: Element_2d.cpp


 ==============================================================================
 void Element2D::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
 {
               KRATOS_TRY
               int FractionalStepNumber = rCurrentProcessInfo[FRACTIONAL_STEP];

               if(FractionalStepNumber == 1) //first step of the fractional step solution 
                  Stage1(rLeftHandSideMatrix,rRightHandSideVector,rCurrentProcessInfo);

               if(FractionalStepNumber == 2) //second step of the fractional step solution 
                  Stage2(rLeftHandSideMatrix,rRightHandSideVector,rCurrentProcessInfo);

               if(FractionalStepNumber == 3) //third step of the fractional step solution 	
                  Stage3(rLeftHandSideMatrix,rRightHandSideVector,rCurrentProcessInfo);

                KRATOS_CATCH("")
  } 
 ==============================================================================
 


In this function we should introduce variable FractionalStepNumber to specify the number of the step that we want to run in this element. The value of this variable will be assigned later with FRACTIONAL_STEP variable in Python.Three functions of Stage1, Stage2 and Stage3 are defined for three steps and each of them contains different system of matrices and equations that we need to solve for that step.

Defining Degrees of Freedom

The next modification is to add the Degrees of Freedom(DOF) for each step because for each step we may have different DOFs or matrix size according to our problem. Suppose that we have following DOFs for each three steps:

Two DOFs for step 1: POLARIZATION_X and POLARIZATION_Y 
One DOF for step 2 : ELECTRIC_POTENTIAL 
Two DOFs for step 3: DISPLACEMENT_X and DISPLACEMENT_Y

then we should add these DOFs to our element. The modifications are in the functions of EquationIdVector and GetDofList in Element_2d.cpp :


Element cpp file: Element_2d.cpp

EquationIdVector function:

 ==============================================================================

   void Element2D::EquationIdVector(EquationIdVectorType& rResult, ProcessInfo& CurrentProcessInfo)
        {
        int number_of_nodes = GetGeometry().size();
        int dim = GetGeometry().WorkingSpaceDimension();
        unsigned int dim2 = number_of_nodes*dim;
        unsigned int FractionalStepNumber = CurrentProcessInfo[FRACTIONAL_STEP];

        if(FractionalStepNumber == 1) //step 1
           {
           if(rResult.size() != dim2)
              rResult.resize(dim2);
              for (int i=0;i<number_of_nodes;i++)
                  {
                  int index = i*dim;
                  rResult[index] = GetGeometry()[i].GetDof(POLARIZATION_X).EquationId();
                  rResult[index + 1] = GetGeometry()[i].GetDof(POLARIZATION_Y).EquationId();
                  }
           }

        if(FractionalStepNumber == 2) //step 2
           {
           if(rResult.size() != number_of_nodes)
              rResult.resize(number_of_nodes,false);	
 	    for (int i=0;i<number_of_nodes;i++)
                rResult[i] = GetGeometry()[i].GetDof(ELECTRIC_POTENTIAL).EquationId();
           }

        if(FractionalStepNumber == 3) //step 3
           {
           if(rResult.size() != dim2)
              rResult.resize(dim2);
           for (int i=0;i<number_of_nodes;i++)
               {
               int index = i*dim;
               rResult[index] = GetGeometry()[i].GetDof(DISPLACEMENT_X).EquationId();
               rResult[index + 1] = GetGeometry()[i].GetDof(DISPLACEMENT_Y).EquationId();
               }
           }
      }

==============================================================================

GetDofList function:

==============================================================================

  void Element2D::GetDofList(DofsVectorType& ElementalDofList,ProcessInfo& CurrentProcessInfo)
       {
       unsigned int FractionalStepNumber = CurrentProcessInfo[FRACTIONAL_STEP];
       ElementalDofList.resize(0);

       if(FractionalStepNumber == 1) //step 1
          {
          for (unsigned int i=0;i<GetGeometry().size();i++)
              {
              ElementalDofList.push_back(GetGeometry()[i].pGetDof(POLARIZATION_X));
              ElementalDofList.push_back(GetGeometry()[i].pGetDof(POLARIZATION_Y));
              }
          }

       if(FractionalStepNumber == 2) //step 2
          {
          for (unsigned int i=0;i<GetGeometry().size();i++)
              ElementalDofList.push_back(GetGeometry()[i].pGetDof(ELECTRIC_POTENTIAL));
          }

       if(FractionalStepNumber == 3) //step 3
          {
          for (unsigned int i=0;i<GetGeometry().size();i++)
              {
              ElementalDofList.push_back(GetGeometry()[i].pGetDof(DISPLACEMENT_X));
              ElementalDofList.push_back(GetGeometry()[i].pGetDof(DISPLACEMENT_Y));
              }
          }
        }

=============================================================================


After adding these DOFs for each step, Kratos automatically referes to them according to FractionalStepNumber, for constructing the RHS and LHS matrices. For the final modification, we should define three step functions which include the definitation of RHS and LHS matrices. These definitions are related to your problem and here we only mention the total structure of these functions. You can define your problem, RHS and LHS matrices according to DOFs for each step similar to the common elements. You also need to define these three functions in the header file (Element_2d.h):

Defining step functions

Element cpp file: Element_2d.cpp

  //The first Equation 
  void Element2D::Stage1(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector,ProcessInfo& rCurrentProcessInfo)
  {
  KRATOS_TRY;
 
  // Define your RHS and LHS matrices for the first problem here 

  KRATOS_CATCH("");
  }


  //The second Equation 
 void ELement2D::Stage2(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector,ProcessInfo& rCurrentProcessInfo)
 {
 KRATOS_TRY;

 // Define your RHS and LHS matrices for the second problem here 

 KRATOS_CATCH("");
 }


 //The third Equation 
 void Element2D::Stage3(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector,ProcessInfo& rCurrentProcessInfo)
 {
 KRATOS_TRY;

 // Define your RHS and LHS matrices for the third problem here 

 KRATOS_CATCH("");
 }


Element header file: Element_2d.h

void Stage1(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo);

void Stage2(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo);

void Stage3(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo);


Python interface files

Now our element files are ready and we should move to Python interface file and add or modify some lines to run our multi-step application. The first thing that we need is to define DOFs in our solver file. We can add all the DOFs ,which we defined for three steps, in AddDofs section of this file which is located in python_scripts folder of the application:

static_Element2D_solver.py (solver file):


def AddDofs(model_part):
  for node in model_part.Nodes:

      #adding dofs
       node.AddDof(POLARIZATION_X);
       node.AddDof(POLARIZATION_Y);
       node.AddDof(ELECTRIC_POTENTIAL);
       node.AddDof(DISPLACEMENT_X);
       node.AddDof(DISPLACEMENT_Y);


The second step is to add the solver to Python interface file in the test_examples folder or any folder that we put Python run codes. The important point is that for each of the mentioned three steps we need separate solver and we should not use the same solver for more than one step because the number of DOFs may be different for two steps. To define three separate solvers we can add following lines to our Python file which we name it Run.py:


Run.py (Run file):


static_Element2D_solver.AddDofs(model_part)
#creating the first solver object for step 1
solver1 = static_Element2D_solver.StaticElement2DSolver(model_part,domain_size)
solver1.Initialize()


#creating the second solver object for step 2
solver2 = static_Element2D_solver.StaticElement2DSolver(model_part,domain_size)
solver2.Initialize()


#creating the third solver object for step 3
solver3 = static_Element2D_solver.StaticElement2DSolver(model_part,domain_size)
solver3.Initialize()


After defining these solvers, you can solve each step of your problem by adding the following lines in your run file:

Run.py (Run file):

#Step 1
#FractionalStepNumber = 1
model_part.ProcessInfo[FRACTIONAL_STEP] = 1
solver1.Solve()
#Step 2
#FractionalStepNumber = 2
model_part.ProcessInfo[FRACTIONAL_STEP] = 2
solver2.Solve()
#Step 3
#FractionalStepNumber = 3
model_part.ProcessInfo[FRACTIONAL_STEP] = 3
solver3.Solve()

where model_part.ProcessInfo[FRACTIONAL_STEP] assigns the current step number. In this example, steps 1 to 3 are consecutively solved but you can change this order according to your problem or put these solver lines in a loop to make a transient or iterative schemes.

Personal tools
Categories