How to compute the residual for Structural Mechanics case ?

[azzeddinetiba]

Description I think my question can be devided in three parts.

1.) How can we retrieve the residual after the solution of a structural mechanics problem from Python. Meaning after running SolveSolutionStep() , can we compute/retrieve the residual after that ?

2.) Can we compute this residual with a prescribed displacement field ? For example, after running InitializeSolutionStep(), and setting the displacement field ourselves, can the residual still be computed the same way ?

3.) Does that hold even for a nonlinear problem with point loads ? (Specifically, I want to apply this to a Hyperlastic problem, with a Total Lagrange Element with point loads) since it's not an essential BC, is there another step needed in that case ?

Thank you in advance !

Scope

• StructuralMechanicsApplication
• Python

[Maziarb]

Hi, it is also my question. I set the displacement in the full order model (FOM) using projected data from the reduced model (ROM). But to compute the residual in FOM I get some difficulties and need help too. Here is the piece of code for computing the residuals:

def compute_residuals(analysis):
    # Access the ComputingModelPart
    model_part = analysis._GetSolver().GetComputingModelPart()
    system_size = model_part.GetCommunicator().GetDataCommunicator().Sum(model_part.NumberOfNodes() * 3, 0)

    residuals = Kratos.Vector(system_size)
    for i in range(system_size):
        residuals[i] = 0.0

    # Get the Scheme from the solver
    scheme = analysis._GetSolver()._GetScheme()

    # Get the BuilderAndSolver instance
    builder_and_solver = analysis._GetSolver()._GetBuilderAndSolver()

    # Assemble the global residual vector
    builder_and_solver.BuildRHS(scheme, model_part, residuals)

    return residuals

Hope someone could give us some solutions to that.

Thanks :)

[NicolasSR]

Hi. I'll answer following @azzeddinetiba's 3-parts structure.

In my case I need the residuals to be saved in a matrix and saved to disk, so I do it in Python. That is the way I will show you. It shouldn't be much different if you want to implement it in C++ and use the residual internally.

1) For the FOM problem, I simply use an analysis stage on top of the application's one, with the following methods:

def ModifyInitialGeometry(self):
    super().ModifyInitialGeometry()

    self.residuals_matrix = []

    self.modelpart = self._GetSolver().GetComputingModelPart()
    self.strategy = self._GetSolver()._GetSolutionStrategy()
    self.buildsol = self._GetSolver()._GetBuilderAndSolver()
    self.scheme = self._GetSolver()._GetScheme()
    self.space = KratosMultiphysics.UblasSparseSpace()

def FinalizeSolutionStep(self):
    super().FinalizeSolutionStep()

    # Initialize residuals vector
    b = self.strategy.GetSystemVector()
    self.space.SetToZeroVector(b)

    # Build residuals based on FOM results, and append to matrix
    self.buildsol.BuildRHS(self.scheme, self.modelpart, b)
    b=np.array(b, copy=False)
    self.residuals_matrix.append(b)

def Finalize(self):
    super().Finalize()

    # Write to disk
    np.save("residuals_matrix.npy", np.array(self.residuals_matrix, copy=False))

This should be pretty much equivalent to @Maziarb's approach.

2) Indeed, you can compute the residuals for a custom displacement as long as you update that displacement before calling the BuildRHS method. The way I do it is changing the FinalizeSolutionStep method above to:

def FinalizeSolutionStep(self):
    super().FinalizeSolutionStep()

    # Initialize residuals vector
    b = self.strategy.GetSystemVector()
    self.space.SetToZeroVector(b)

    # Update displacements and pointlads to our own
    self.project_displacements(my_displacements, my_pointloads)

    # Build residuals based on our displacements, and append to matrix
    self.buildsol.BuildRHS(self.scheme, self.modelpart, b)
    b=np.array(b, copy=False)
    self.residuals_matrix.append(b)

def project_displacements(self, my_displacements, my_pointloads):

    dim = 2
    nodes_array=self.modelpart.Nodes

    # Update positions and displacement variables
    x0_vec = self.var_utils.GetInitialPositionsVector(nodes_array,dim)
    self.var_utils.SetSolutionStepValuesVector(nodes_array, KratosMultiphysics.DISPLACEMENT, my_displacements, 0)
    x_vec=x0_vec+my_displacements
    self.var_utils.SetCurrentPositionsVector(nodes_array,x_vec)

    # Update point loads
    conditions_array=self.modelpart.Conditions
    for i, condition in enumerate(conditions_array):
        condition.SetValue(StructuralMechanicsApplication.POINT_LOAD, my_pointloads[i])

3) As can be seen in the code above, if we apply point or line loads, we must update the corresponding conditions accordingly. Otherwise those forces will appear added on the residual. I use this approach on non-linear problems and it holds.

Hope this helped! And don't hesitate to ask further

[azzeddinetiba]

Thank you @NicolasSR . Unfortunately this doesn't seem to work with Point Loads. (or At least as they are configured in my case). So with just computing the residuals after the StructuralMechancis solution, the modelpart conditions seem to have 0 loads even if they are prescribed at POINT_LOAD . But it does give a very small residual. However when prescribing my own displacement (even an accurate one) the residuals are large.

So, for the StructuralMechanics solution, this :

  def SolveSolutionStep(self):

      # ========= To compute the residuals ===================
      model_part = self._analysis_stage._GetSolver().GetComputingModelPart()
      system_size = model_part.GetCommunicator().GetDataCommunicator().Sum(model_part.NumberOfNodes() * 2, 0)
      residuals = KM.Vector(system_size)
      for i in range(system_size):
          residuals[i] = 0.0

      scheme = self._analysis_stage._GetSolver()._GetScheme()

      builder_and_solver = self._analysis_stage._GetSolver()._GetBuilderAndSolver()
      conditions_array=model_part.Conditions

      print("\n Prescribed LOAD norm", np.linalg.norm(KM.VariableUtils().GetSolutionStepValuesVector(
                      self._analysis_stage._GetSolver().GetComputingModelPart().Nodes, 
                      StructuralMechanicsApplication.POINT_LOAD, 0, 2)), " \n")
      for i, condition in enumerate(conditions_array):
          thisId = condition.GetNode(0).Id
          print("\n cnd PT LOAD ", thisId, " ", condition.GetValue(StructuralMechanicsApplication.POINT_LOAD), "\n")

      # ======= Actual solution ============
      super().SolveSolutionStep()

      builder_and_solver.BuildRHS(scheme, model_part, residuals)
      print("\n", np.linalg.norm(np.asarray(residuals)), "\n")

gives

 Prescribed LOAD norm 111.72215732198302  

 cnd PT LOAD  39   [3](0,0,0) 

 cnd PT LOAD  41   [3](0,0,0) 

 cnd PT LOAD  46   [3](0,0,0) 

 cnd PT LOAD  49   [3](0,0,0) 

 cnd PT LOAD  53   [3](0,0,0) 
...

 7.721769862193473e-06 

Is it normal that the conditions are empty ?

If I prescribe a given displacement, i.e

#super().SolveSolutionStep()              
KM.VariableUtils().SetSolutionStepValuesVector(self._analysis_stage._GetSolver().GetComputingModelPart().Nodes,
                                                               KM.DISPLACEMENT, myDisplacement, 0)
x0_vec = KM.VariableUtils().GetInitialPositionsVector(self._analysis_stage._GetSolver().GetComputingModelPart().Nodes,2)
x_vec=x0_vec+predicted_disp
KM.VariableUtils().SetCurrentPositionsVector(self._analysis_stage._GetSolver().GetComputingModelPart().Nodes,x_vec)

It gives something like

 Prescribed LOAD norm 111.60857092935072  

 cnd PT LOAD  39   [3](0,0,0) 

 cnd PT LOAD  41   [3](0,0,0) 

 cnd PT LOAD  46   [3](0,0,0) 

 cnd PT LOAD  49   [3](0,0,0) 
...

 20.61074378468338 

Do you guys have any idea ? PS: This is done inside a CoSimulationApplication, the above functions are made on top of those of StructuralMechanicsWrapper . Is something happening on the InitializeSolutionStep() ?

Thank you !

[NicolasSR]

Hello again,

Are you sure the point loads are being applied correctly? Can you check if the displacements of your model are actually being affected by the load?

In my case I get something like:

Prescribed LOAD norm 0.0  

cnd PT LOAD  10   [3](0,-600,0) 

2.9351115032693715e-05

My guess is you're trying to apply the point load as a node variable, instead of a condition, and it isn't actually taken into account to solve the step. That would explain why the first norm is higher than zero for you, while the values of your conditions are all zeros.

When prescribing your own displacement, it will no longer match the point loads being applied to your model (unless you modify the loads yourself with the appropriate values) and therefore you will have high values for the residual at the indexes corresponding to the nodes with the point loads. That explains the residual norm being 20 in your second case.

[azzeddinetiba]

Hi,

The displacement field does get updated each time and has the expected results, according to the prescribed loads. I can retrieve it via GetSolutionStepValuesVector()

So, - if I understood you correctly- this doesn't explain why the residuals are 0 when solved from Kratos and only non-zero when I prescribe my own displacement.

The point loads are actually applied through the Map() function of the MappingApplication

[RiccardoRossi]

my guess is that you are setting the POINT_LOAD in the historical database instead of on the nonhistorical one (or the contrary)

[azzeddinetiba]

Hello @RiccardoRossi , Thank you for replying,

I'm using the MappingApplication 'Nearest neighbors' , and as used in the tutorials, the default location "historical" is used https://github.com/KratosMultiphysics/Kratos/blob/f8aeb637293ce5177d6bca41d629fd31a8b4a24b/applications/CoSimulationApplication/python_scripts/coupling_interface_data.py#L67

But I'm not sure how (and at which point in the code) the mapper assigns the values of the variables, and if it does indeed on the historical or the nonhistorical database.

But If I understood correctly, this should affect the Kratos results as well and not just when I prescribe an external displacement field ?