Implementing thermal solver

The following tutorial will explain how to implement a Solver from scratch, in this case applied to the particular case of the thermal problem, following the lines of the already presented on the other tutorials. We will skip the most advanced points of the construction of the Solver and we will create a basic working solver for the sake of academic purposes, so the resulting file will not coincide with the real solver on the repository.

Imports

The solver as any element of Kratos requires to import the corresponding libraries and applications. Like out objective on mind is to create a thermal problem, we will import the base KratosMultiphysics library as well as the ConvectionDiffusionApplication. For importing files from the filesystem we will import os python library, which will be helpful.

from __future__ import print_function, absolute_import, division  # makes KratosMultiphysics backward compatible with python 2.6 and 2.7

# Importing the Kratos Library
import KratosMultiphysics

# Check that applications were imported in the main script
KratosMultiphysics.CheckRegisteredApplications("ConvectionDiffusionApplication")

# Import applications
import KratosMultiphysics.ConvectionDiffusionApplication as ConvectionDiffusionApplication

# Other imports
import os

Constructing the solver

We will create a Solver, following was is done in other applications Solver, for the sake of consistency, so in first place we define the function CreateSolver, which is common among all the solvers and therefore it is nececessary to be called that way. This function will use as input a ModelPart and confuguration Parameters:

def CreateSolver(main_model_part, custom_settings):
    return ConvectionDiffusionSolver(main_model_part, custom_settings)

As we see this function call a class called ConvectionDiffusionSolver, so we need to define our solver in first place, for that we define the constructor or __init__ function:

class ConvectionDiffusionSolver(object):
    """The base class for convection-diffusion solvers.
    This class provides functions for importing and exporting models,
    adding nodal variables and dofs and solving each solution step.
    """
    def __init__(self, main_model_part, custom_settings):
        default_settings = KratosMultiphysics.Parameters("""
        {
            "echo_level": 0,
            "buffer_size": 2,
            "model_import_settings": {
                "input_type": "mdpa",
                "input_filename": "unknown_name"
            },
            "computing_model_part_name" : "Thermal",
            "material_import_settings" :{
                "materials_filename": ""
            },
            "clear_storage": false,
            "residual_relative_tolerance": 1.0e-4,
            "residual_absolute_tolerance": 1.0e-9,
            "max_iteration": 10,
            "linear_solver_settings":{
                "solver_type": "BICGSTABSolver",
                "preconditioner_type": "DiagonalPreconditioner",
                "max_iteration": 5000,
                "tolerance": 1e-9,
                "scaling": false
            },
            "element_replace_settings" : {
                "element_name" : "EulerianConvDiff",
                "condition_name" : "Condition"
            },
            "problem_domain_sub_model_part_list": ["conv_diff_body"],
            "processes_sub_model_part_list": [""]
        }
        """)
    
        # Overwrite the default settings with user-provided parameters.
        self.settings = custom_settings
        self.settings.ValidateAndAssignDefaults(default_settings)

        self.main_model_part = main_model_part
        KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionSolver]:: ", "Construction finished")

First we define the default Parameters, which will overwrite the not define configurations of our custom_settings, this is done with the ValidateAndAssignDefaults method. We define the main_model_part as a part of the solver, using the self to define it as an instance attribute.

Common methods for the solver

The following are common methods defined in all solvers in order to have an interoperability between them to couple different physical problems.

Add variables

The AddVariables is the responsible on adding the nodal historical variables to the model part. We will need to define a buffer in order to be able to access this historical variables. We will see how in following sections.

We add the variables corresponding to the standard thermal convection-diffusion problem. We will need to add this variables to the CONVECTION_DIFFUSION_SETTINGS (this is something specific of the ConvectionDiffusionAplication, and is done in order to solve different types convection diffusion problems).

def AddVariables(self):
    ''' Add nodal solution step variables ( later add to CONVECTION_DIFFUSION_SETTINGS)
    '''
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DENSITY)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.CONDUCTIVITY)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.TEMPERATURE)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.HEAT_FLUX)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.FACE_HEAT_FLUX)
    self.main_model_part.AddNodalSolutionStepVariable(ConvectionDiffusionApplication.PROJECTED_SCALAR1)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.CONVECTION_VELOCITY)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.MESH_VELOCITY)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.VELOCITY)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.SPECIFIC_HEAT)
    self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.REACTION_FLUX)

    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]:: ", "Variables ADDED")

Add DoFs

In this case we add to the nodes the degrees of freedom corresponding to the problems we want to solve, in our case only the TEMPERATURE.

def AddDofs(self):
        KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.TEMPERATURE, KratosMultiphysics.REACTION_FLUX, self.main_model_part)
        KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]:: ", "DOF's ADDED")

We the VariableUtils class, which allows us to do this operation in parallel, in case we want to use the methods available directly on the Node class we can do:

def AddDofs(self):
        for node in self.main_model_part.Nodes:
            node.AddDof(KratosMultiphysics.TEMPERATURE, KratosMultiphysics.REACTION_FLUX)  
        KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]:: ", "DOF's ADDED")  

Read external file (*.mdpa file)

def ReadModelPart(self):
    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]::", "Reading model part.")
    problem_path = os.getcwd()
    input_filename = self.settings["model_import_settings"]["input_filename"].GetString()
    # Import model part from mdpa file.
    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]::", "Reading model part from file: " + os.path.join(problem_path, input_filename) + ".mdpa")
    KratosMultiphysics.ModelPartIO(input_filename).ReadModelPart(self.main_model_part)
    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]::", "Finished reading model part from mdpa file.")
    KratosMultiphysics.Logger.PrintInfo("ModelPart", self.main_model_part)
    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]:: ", "Finished reading model part.")

Prepare model part

    def _execute_after_reading(self):
        """Prepare computing model part and import constitutive laws. """
        # Auxiliary parameters object for the CheckAndPepareModelProcess
        params = KratosMultiphysics.Parameters("{}")
        params.AddValue("computing_model_part_name",self.settings["computing_model_part_name"])
        params.AddValue("problem_domain_sub_model_part_list",self.settings["problem_domain_sub_model_part_list"])
        params.AddValue("processes_sub_model_part_list",self.settings["processes_sub_model_part_list"])
        # Assign mesh entities from domain and process sub model parts to the computing model part.
        import check_and_prepare_model_process_convection_diffusion as check_and_prepare_model_process
        check_and_prepare_model_process.CheckAndPrepareModelProcess(self.main_model_part, params).Execute()

        # Import constitutive laws.
        materials_imported = self.import_materials()
        if materials_imported:
            self.print_on_rank_zero("::[ConvectionDiffusionBaseSolver]:: ", "Materials were successfully imported.")
        else:
            self.print_on_rank_zero("::[ConvectionDiffusionBaseSolver]:: ", "Materials were not imported.")

    def _set_and_fill_buffer(self):
        """Prepare nodal solution step data containers and time step information. """
        # Set the buffer size for the nodal solution steps data. Existing nodal
        # solution step data may be lost.
        required_buffer_size = self.settings["buffer_size"].GetInt()
        if required_buffer_size < self.GetMinimumBufferSize():
            required_buffer_size = self.GetMinimumBufferSize()
        current_buffer_size = self.main_model_part.GetBufferSize()
        buffer_size = max(current_buffer_size, required_buffer_size)
        self.main_model_part.SetBufferSize(buffer_size)
        # Cycle the buffer. This sets all historical nodal solution step data to
        # the current value and initializes the time stepping in the process info.
        delta_time = self.main_model_part.ProcessInfo[KratosMultiphysics.DELTA_TIME]
        time = self.main_model_part.ProcessInfo[KratosMultiphysics.TIME]
        step =-buffer_size
        time = time - delta_time * buffer_size
        self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.TIME, time)
        for i in range(0, buffer_size):
            step = step + 1
            time = time + delta_time
            self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.STEP, step)
            self.main_model_part.CloneTimeStep(time)

    def _get_element_condition_replace_settings(self):
        num_nodes_elements = 0
        if (len(self.main_model_part.Elements) > 0):
            num_nodes_elements = len(self.main_model_part.Elements[1].GetNodes())

        ## Elements
        if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2:
            if (num_nodes_elements == 3):
                self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff2D")
            else:
                self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff2D4N")
        elif self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3:
            if (self.settings["element_replace_settings"]["element_name"].GetString() == "EulerianConvDiff"):
                self.settings["element_replace_settings"]["element_name"].SetString("EulerianConvDiff3D")
        else:
            raise Exception("DOMAIN_SIZE not set")
        
        ## Conditions
        num_nodes_conditions = 0
        if (len(self.main_model_part.Conditions) > 0):
            num_nodes_conditions = len(self.main_model_part.Conditions[1].GetNodes())
        if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2:
            self.settings["element_replace_settings"]["condition_name"].SetString("LineCondition2D2N")
        elif self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3:
            self.settings["element_replace_settings"]["condition_name"].SetString("SurfaceCondition3D3N")
        else:
            raise Exception("DOMAIN_SIZE not set")

        return self.settings["element_replace_settings"]

Initialize method

def Initialize(self):
    """Perform initialization after adding nodal variables and dofs to the main model part. """
    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]:: ", "Initializing ...")

    # We define the convection diffusion settings
    convention_diffusion_settings = KratosMultiphysics.ConvectionDiffusionSettings()
    convention_diffusion_settings.SetDensityVariable(KratosMultiphysics.DENSITY)
    convention_diffusion_settings.SetDiffusionVariable(KratosMultiphysics.CONDUCTIVITY)
    convention_diffusion_settings.SetUnknownVariable(KratosMultiphysics.TEMPERATURE)
    convention_diffusion_settings.SetVolumeSourceVariable(KratosMultiphysics.HEAT_FLUX)
    convention_diffusion_settings.SetSurfaceSourceVariable(KratosMultiphysics.FACE_HEAT_FLUX)
    convention_diffusion_settings.SetProjectionVariable(ConvectionDiffusionApplication.PROJECTED_SCALAR1)
    convention_diffusion_settings.SetConvectionVariable(KratosMultiphysics.CONVECTION_VELOCITY)
    convention_diffusion_settings.SetMeshVelocityVariable(KratosMultiphysics.MESH_VELOCITY)
    convention_diffusion_settings.SetVelocityVariable(KratosMultiphysics.VELOCITY)
    convention_diffusion_settings.SetSpecificHeatVariable(KratosMultiphysics.SPECIFIC_HEAT)
    convention_diffusion_settings.SetReactionVariable(KratosMultiphysics.REACTION_FLUX)
        
    self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.CONVECTION_DIFFUSION_SETTINGS, convention_diffusion_settings)
        
    # The mechanical solution strategy is created here if it does not already exist.
    computing_model_part = self.main_model_part.GetSubModelPart(self.settings["computing_model_part_name"].GetString())
    conv_diff_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticScheme()
    import linear_solver_factory
    linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"])
    conv_diff_convergence_criterion = KratosMultiphysics.ResidualCriteria(self.settings["residual_relative_tolerance"].GetDouble(), self.settings["residual_absolute_tolerance"].GetDouble())
    builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(linear_solver)
    self.conv_diff_strategy = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(computing_model_part,  conv_diff_scheme, linear_solver, conv_diff_convergence_criterion, builder_and_solver, self.settings["max_iteration"].GetInt(), True, False, False)
    self.conv_diff_strategy.SetEchoLevel(self.settings["echo_level"].GetInt())
    self.conv_diff_strategy.Initialize()
    self.conv_diff_strategy.Check()
    if self.settings["clear_storage"].GetBool():
        self.conv_diff_strategy.Clear()
    KratosMultiphysics.Logger.PrintInfo("::[ConvectionDiffusionBaseSolver]:: ", "Finished initialization.")

Solving methods

The Solver as well as the strategy, has the same solving steps as the strategies (Newton-Rapshon and so on), this methods are InitializeSolutionStep, Predict, SolveSolutionStep and FinalizeSolutionStep, which can be call all together in only one method Solve.

    def Solve(self):
        if self.settings["clear_storage"].GetBool():
            self.conv_diff_strategy.Clear()
        self.conv_diff_strategy.Solve()

    def InitializeSolutionStep(self):
        self.conv_diff_strategy.InitializeSolutionStep()

    def Predict(self):
        self.conv_diff_strategy.Predict()

    def SolveSolutionStep(self):
        is_converged = self.conv_diff_strategy.SolveSolutionStep()
        return is_converged

    def FinalizeSolutionStep(self):
        self.conv_diff_strategy.FinalizeSolutionStep()

Reading materials

    def import_materials(self):
        materials_filename = self.settings["material_import_settings"]["materials_filename"].GetString()
        if (materials_filename != ""):
            import read_materials_process
            # Create a dictionary of model parts.
            Model = KratosMultiphysics.Model()
            Model.AddModelPart(self.main_model_part)
            # Add constitutive laws and material properties from json file to model parts.
            read_materials_process.ReadMaterialsProcess(Model, self.settings["material_import_settings"])
            
            # We set the properties that are nodal
            self._assign_nodally_properties()
            
            materials_imported = True
        else:
            materials_imported = False
        return materials_imported

    def _assign_nodally_properties(self):
        
        # We transfer the values of the con.diff variables to the nodes
        with open(self.settings["material_import_settings"]["materials_filename"].GetString(), 'r') as parameter_file:
            materials = KratosMultiphysics.Parameters(parameter_file.read())
            
        for i in range(materials["properties"].size()):
            model_part = self.main_model_part.GetSubModelPart(materials["properties"][i]["model_part_name"].GetString())
            mat = materials["properties"][i]["Material"]
            
            for key, value in mat["Variables"].items():
                var = KratosMultiphysics.KratosGlobals.GetVariable(key)
                if (self._check_variable_to_set(var)):
                    if value.IsDouble():
                        KratosMultiphysics.VariableUtils().SetScalarVar(var, value.GetDouble(), model_part.Nodes)
                    elif value.IsVector():
                        KratosMultiphysics.VariableUtils().SetVectorVar(var, value.GetVector(), model_part.Nodes)
                    else:
                        raise ValueError("Type of value is not available")