Thermally coupled Constitutive Laws

[AlejandroCornejo]

Description Dear @KratosMultiphysics/structural-mechanics and @KratosMultiphysics/constitutive-models

I'm implementing all the code infraestructure in order to be able to simulate thermo-mechanical calculations with a solver that couples the Structural and the ConvectionDiffusion Apps.

The solver was already there long time ago, I've just improved it a bit and corrected a detail.

This is a work in progress PR and I'd like to reach consensus during its implementation.

My main concerns are:

1. In order to reduce code repetition in the thermal_constitutive_laws I'm going to slightly modify the original CL. This change consists in subtituting the rValues.GetProperties()[YOUNG_MODULUS] by a new virtual method called GetMaterialProperty(YOUNG_MODULUS), which in the original cl will do the same operation as before. In the case of the thermal laws, the only method to be modified would be the GetMaterialProperty(rVariable, rValues), where we modify the property according to a certain TEMPERATURE Table.
2. The setting of the REFERENCE_TEMPERATURE, different options are implemented: using the GetValue() in the Geometry (has to be set outside by someone) and being a material property as:

void ThermalLinearPlaneStrain::InitializeMaterial(
    const Properties& rMaterialProperties,
    const GeometryType& rElementGeometry,
    const Vector& rShapeFunctionsValues
    )
{
    if (rElementGeometry.Has(REFERENCE_TEMPERATURE)) {
        mReferenceTemperature = rElementGeometry.GetValue(REFERENCE_TEMPERATURE);
    } else if (rMaterialProperties.Has(REFERENCE_TEMPERATURE)) {
        mReferenceTemperature = rMaterialProperties[REFERENCE_TEMPERATURE];
    }
}

Changelog Please summarize the changes in one list to generate the changelog: E.g.

• Erased the convective term in the Eulerian calculation
• Adding a virtual method to the original solid laws for retrieving material properties
• Retrieving material properties in the original CL by means of the GetMaterialProperty(rVariable, rValues) method
• Creation of the new thermal CL inside the ConstitutiveLawsApp, only the elastic basic ones
• Creation of an empty utility container of the ones to be exported to Python
• etc...

Documentation

Now, the materials.json for the structural problem looks as this:

{
    "properties" : [{
        "model_part_name" : "Structure.Parts_Solid_Solid_Auto1",
        "properties_id"   : 1,
        "Material"        : {
            "constitutive_law" : {
                "name" : "ThermalLinearPlaneStrain"
                
            },
            "Variables"        : {
                "DENSITY"       : 7850.0,
                "YOUNG_MODULUS" : 2.1e11,
                "POISSON_RATIO" : 0.29,
                "THERMAL_EXPANSION_COEFFICIENT" : 7.2e-6,
                "REFERENCE_TEMPERATURE" : 22.5
            },
            "Tables"           : {
                "TEMPERATURE_vs_E" : {
                        "input_variable"  : "TEMPERATURE",
                        "output_variable" : "YOUNG_MODULUS",
                        "data"            : [[-1,   2.1e11],
                                             [3500, 2.1e8 ],
                                             [1e6,  2.1e8 ]]
                    }
            }
        }
    }]
}

[loumalouomega]

1. In order to reduce code repetition in the `thermal_constitutive_laws `I'm going to slightly modify the original CL. This change consists in subtituting the `rValues.GetProperties()[YOUNG_MODULUS]` by a new virtual method called `GetMaterialProperty(YOUNG_MODULUS)`, which in the original cl will do the same operation as before. In the case of the thermal laws, the only method to be modified would be the `GetMaterialProperty(rVariable)`, where we modifiy the property according to a certain `TEMPERATURE` Table.

This has a small overhead..., I am OK because there is not easy manner. Except maybe templated static auxiliary methods...

[AlejandroCornejo]

1. In order to reduce code repetition in the `thermal_constitutive_laws `I'm going to slightly modify the original CL. This change consists in subtituting the `rValues.GetProperties()[YOUNG_MODULUS]` by a new virtual method called `GetMaterialProperty(YOUNG_MODULUS)`, which in the original cl will do the same operation as before. In the case of the thermal laws, the only method to be modified would be the `GetMaterialProperty(rVariable)`, where we modifiy the property according to a certain `TEMPERATURE` Table.

This has a small overhead..., I am OK because there is not easy manner. Except maybe templated static auxiliary methods...

Te other option is to repeat all the code and leave alone the original CL... IDK if more templates.. hehe

[AlejandroCornejo]

can you please make some performance tests for adding this virtual method? This is called a lot of times

Other than that the changes look ok to me, but see my question about the access of the props from the elements (which was one of the main reasons #2349 failed)

Sure, before merging we will asses the performance of the implementation

[loumalouomega]

@AlejandroCornejo a minimal C++ test with a minimal hexa could be interesting

[AlejandroCornejo]

@AlejandroCornejo a minimal C++ test with a minimal hexa could be interesting

Yes! a set of test is coming :D

[AlejandroCornejo]

@AlejandroCornejo a minimal C++ test with a minimal hexa could be interesting

I've added several test, one per each added thermal cl in which the TEMPERATURE is involved and tables interpolation of material properties as well.

[AlejandroCornejo]

@philbucher I've been doing the performance test of the original LinearElastic3DLaw with or without the modification.

The test consists in calling 10,000,000 times the CalculateMaterialResponseCauchy method with the following results:

With the virtual method: 5.5625 s Without the virtual metho: 5.4375 s.

[philbucher]

@philbucher I've been doing the performance test of the original LinearElastic3DLaw with or without the modification.

The test consists in calling 10,000,000 times the CalculateMaterialResponseCauchy method with the following results:

With the if: 5.5625 s Without the if: 5.4375 s.

cool can you share the test code? I would like to see it to make sure the compiler didn't optimize anything away

[AlejandroCornejo]

@philbucher I've been doing the performance test of the original LinearElastic3DLaw with or without the modification. The test consists in calling 10,000,000 times the CalculateMaterialResponseCauchy method with the following results: With the if: 5.5625 s Without the if: 5.4375 s.

cool can you share the test code? I would like to see it to make sure the compiler didn't optimize anything away

Sure!!

The code used is:

import KratosMultiphysics as KM
import KratosMultiphysics.StructuralMechanicsApplication as SMA
import KratosMultiphysics.ConstitutiveLawsApplication as CLA
import applications.StructuralMechanicsApplication.tests.test_constitutive_law as test
from time import process_time

def CreateGeometry(model_part, dim):
    # Create new nodes
    node1 = model_part.CreateNewNode(1, 0.0, 0.0, 0.0)
    node2 = model_part.CreateNewNode(2, 1.0, 0.0, 0.0)
    node3 = model_part.CreateNewNode(3, 0.0, 1.0, 0.0)

    if (dim == 2):
        nnodes = 3

        # Allocate a geometry
        geom = KM.Triangle2D3(node1,node2,node3)
    elif (dim == 3):
        nnodes = 4
        node4 = model_part.CreateNewNode(4, 0.0, 0.0, 1.0)

        # Allocate a geometry
        geom = KM.Tetrahedra3D4(node1,node2,node3,node4)
    else:
        raise Exception("Error: bad dimension value: ", dim)
    return [geom, nnodes]

def _set_cl_parameters(cl_options, F, detF, strain_vector, stress_vector,
                       constitutive_matrix, N, DN_DX, model_part, properties, geom):
    # Setting the parameters - note that a constitutive law may not need them all!
    cl_params = KM.ConstitutiveLawParameters()
    cl_params.SetOptions(cl_options)
    cl_params.SetDeformationGradientF(F)
    cl_params.SetDeterminantF(detF)
    cl_params.SetStrainVector(strain_vector)
    cl_params.SetStressVector(stress_vector)
    cl_params.SetConstitutiveMatrix(constitutive_matrix)
    cl_params.SetShapeFunctionsValues(N)
    cl_params.SetShapeFunctionsDerivatives(DN_DX)
    cl_params.SetProcessInfo(model_part.ProcessInfo)
    cl_params.SetMaterialProperties(properties)
    cl_params.SetElementGeometry(geom)

    ## Do all sort of checks
    cl_params.CheckAllParameters() # Can not use this until the geometry is correctly exported to python
    cl_params.CheckMechanicalVariables()
    cl_params.CheckShapeFunctions()
    return cl_params

def CalculateStressFromStrain(StrainVector):
    # cl = CLA.SmallStrainIsotropicPlasticity3DVonMisesVonMises()
    cl = SMA.LinearElastic3DLaw()

    current_model = KM.Model()
    mdpa = current_model.CreateModelPart("test")

    # properties = mdpa.CreateProperties(1)
    properties = KM.Properties(1)

    properties.SetValue(KM.YIELD_STRESS,275E6)
    properties.SetValue(KM.FRACTURE_ENERGY,1e8)
    properties.SetValue(CLA.HARDENING_CURVE,0)
    properties.SetValue(KM.YOUNG_MODULUS,210e9)
    properties.SetValue(KM.POISSON_RATIO,0.3)
    properties.SetValue(KM.CONSTITUTIVE_LAW, cl)

    [geom, nnodes] = CreateGeometry(mdpa, 3)

    N = KM.Vector(nnodes)
    DN_DX = KM.Matrix(nnodes, 3)

    cl.Check(properties, geom, mdpa.ProcessInfo)

    cl_options = KM.Flags()
    cl_options.Set(KM.ConstitutiveLaw.USE_ELEMENT_PROVIDED_STRAIN, True)
    cl_options.Set(KM.ConstitutiveLaw.COMPUTE_CONSTITUTIVE_TENSOR, True)
    cl_options.Set(KM.ConstitutiveLaw.COMPUTE_STRESS, True)

    # Define deformation gradient
    F = KM.Matrix(cl.GetStrainSize(), cl.GetStrainSize())
    detF = 1.0

    stress_vector = KM.Vector(cl.GetStrainSize())
    strain_vector = KM.Vector(cl.GetStrainSize())
    constitutive_matrix = KM.Matrix(cl.GetStrainSize(),cl.GetStrainSize())
    for i in range(0, cl.GetStrainSize()):
        stress_vector[i] = 0.0
        strain_vector[i] = StrainVector[i]
        for j in range(0, cl.GetStrainSize()):
            constitutive_matrix[i,j] = 0.0

    # Setting the parameters - note that a constitutive law may not need them all!
    cl_params = _set_cl_parameters(cl_options, F, detF, strain_vector, stress_vector, constitutive_matrix, N, DN_DX, mdpa, properties, geom)
    cl.InitializeMaterial(properties, geom, N)

    cl.InitializeMaterial(properties, geom, KM.Vector(cl.GetStrainSize()))

    t1_start = process_time()
    times = 10000000
    for i in range(0, times):
        cl.CalculateMaterialResponseCauchy(cl_params)
    t1_stop = process_time()
    print("Elapsed time for performing the CalculateMaterialResponseCauchy " + str(times) + " times " + "in seconds:",t1_stop-t1_start)

    cl.FinalizeMaterialResponseCauchy(cl_params)
    return cl_params, cl



###########################
strain_vector = [0.001, 0.01, 0.001, 0.0, 0.0, 0.0]
values, cl  = CalculateStressFromStrain(strain_vector)

[philbucher]

looks good to me :+1: @RiccardoRossi you also ok?

[loumalouomega]

@philbucher I've been doing the performance test of the original LinearElastic3DLaw with or without the modification. The test consists in calling 10,000,000 times the CalculateMaterialResponseCauchy method with the following results: With the if: 5.5625 s Without the if: 5.4375 s.

cool can you share the test code? I would like to see it to make sure the compiler didn't optimize anything away

Sure!!

The code used is:

import KratosMultiphysics as KM
import KratosMultiphysics.StructuralMechanicsApplication as SMA
import KratosMultiphysics.ConstitutiveLawsApplication as CLA
import applications.StructuralMechanicsApplication.tests.test_constitutive_law as test
from time import process_time

def CreateGeometry(model_part, dim):
    # Create new nodes
    node1 = model_part.CreateNewNode(1, 0.0, 0.0, 0.0)
    node2 = model_part.CreateNewNode(2, 1.0, 0.0, 0.0)
    node3 = model_part.CreateNewNode(3, 0.0, 1.0, 0.0)

    if (dim == 2):
        nnodes = 3

        # Allocate a geometry
        geom = KM.Triangle2D3(node1,node2,node3)
    elif (dim == 3):
        nnodes = 4
        node4 = model_part.CreateNewNode(4, 0.0, 0.0, 1.0)

        # Allocate a geometry
        geom = KM.Tetrahedra3D4(node1,node2,node3,node4)
    else:
        raise Exception("Error: bad dimension value: ", dim)
    return [geom, nnodes]

def _set_cl_parameters(cl_options, F, detF, strain_vector, stress_vector,
                       constitutive_matrix, N, DN_DX, model_part, properties, geom):
    # Setting the parameters - note that a constitutive law may not need them all!
    cl_params = KM.ConstitutiveLawParameters()
    cl_params.SetOptions(cl_options)
    cl_params.SetDeformationGradientF(F)
    cl_params.SetDeterminantF(detF)
    cl_params.SetStrainVector(strain_vector)
    cl_params.SetStressVector(stress_vector)
    cl_params.SetConstitutiveMatrix(constitutive_matrix)
    cl_params.SetShapeFunctionsValues(N)
    cl_params.SetShapeFunctionsDerivatives(DN_DX)
    cl_params.SetProcessInfo(model_part.ProcessInfo)
    cl_params.SetMaterialProperties(properties)
    cl_params.SetElementGeometry(geom)

    ## Do all sort of checks
    cl_params.CheckAllParameters() # Can not use this until the geometry is correctly exported to python
    cl_params.CheckMechanicalVariables()
    cl_params.CheckShapeFunctions()
    return cl_params

def CalculateStressFromStrain(StrainVector):
    # cl = CLA.SmallStrainIsotropicPlasticity3DVonMisesVonMises()
    cl = SMA.LinearElastic3DLaw()

    current_model = KM.Model()
    mdpa = current_model.CreateModelPart("test")

    # properties = mdpa.CreateProperties(1)
    properties = KM.Properties(1)

    properties.SetValue(KM.YIELD_STRESS,275E6)
    properties.SetValue(KM.FRACTURE_ENERGY,1e8)
    properties.SetValue(CLA.HARDENING_CURVE,0)
    properties.SetValue(KM.YOUNG_MODULUS,210e9)
    properties.SetValue(KM.POISSON_RATIO,0.3)
    properties.SetValue(KM.CONSTITUTIVE_LAW, cl)

    [geom, nnodes] = CreateGeometry(mdpa, 3)

    N = KM.Vector(nnodes)
    DN_DX = KM.Matrix(nnodes, 3)

    cl.Check(properties, geom, mdpa.ProcessInfo)

    cl_options = KM.Flags()
    cl_options.Set(KM.ConstitutiveLaw.USE_ELEMENT_PROVIDED_STRAIN, True)
    cl_options.Set(KM.ConstitutiveLaw.COMPUTE_CONSTITUTIVE_TENSOR, True)
    cl_options.Set(KM.ConstitutiveLaw.COMPUTE_STRESS, True)

    # Define deformation gradient
    F = KM.Matrix(cl.GetStrainSize(), cl.GetStrainSize())
    detF = 1.0

    stress_vector = KM.Vector(cl.GetStrainSize())
    strain_vector = KM.Vector(cl.GetStrainSize())
    constitutive_matrix = KM.Matrix(cl.GetStrainSize(),cl.GetStrainSize())
    for i in range(0, cl.GetStrainSize()):
        stress_vector[i] = 0.0
        strain_vector[i] = StrainVector[i]
        for j in range(0, cl.GetStrainSize()):
            constitutive_matrix[i,j] = 0.0

    # Setting the parameters - note that a constitutive law may not need them all!
    cl_params = _set_cl_parameters(cl_options, F, detF, strain_vector, stress_vector, constitutive_matrix, N, DN_DX, mdpa, properties, geom)
    cl.InitializeMaterial(properties, geom, N)

    cl.InitializeMaterial(properties, geom, KM.Vector(cl.GetStrainSize()))

    t1_start = process_time()
    times = 10000000
    for i in range(0, times):
        cl.CalculateMaterialResponseCauchy(cl_params)
    t1_stop = process_time()
    print("Elapsed time for performing the CalculateMaterialResponseCauchy " + str(times) + " times " + "in seconds:",t1_stop-t1_start)

    cl.FinalizeMaterialResponseCauchy(cl_params)
    return cl_params, cl



###########################
strain_vector = [0.001, 0.01, 0.001, 0.0, 0.0, 0.0]
values, cl  = CalculateStressFromStrain(strain_vector)

With C++17 we can use if constexpr...

[philbucher]

hm now I am again a bit lost

@AlejandroCornejo can you please point again to the place in the code with the if? Thx

[AlejandroCornejo]

okay! As far as I know both @loumalouomega and @philbucher agree with the current state of the PR. What about @RiccardoRossi and @rubenzorrilla ? Regards

[rubenzorrilla]

I've not much to add in here. But in any case, this is only adding new features that do not interfere with the existent ones (we checked that the virtual impact is negligible). Hence, I'd go ahead.

[adityaghantasala]

A general question : Are you planning to add an example case to the Examples repository ?

[AlejandroCornejo]

A general question : Are you planning to add an example case to the Examples repository ?

I can definitely do it!

[adityaghantasala]

A general question : Are you planning to add an example case to the Examples repository ?

I can definitely do it!

Thank you that will be really helpful.

[AlejandroCornejo]

A general question : Are you planning to add an example case to the Examples repository ?

I can definitely do it!

Thank you that will be really helpful.

added in https://github.com/KratosMultiphysics/Examples/pull/88

[AlejandroCornejo]

ping @RiccardoRossi

[AlejandroCornejo]

After talking with @RiccardoRossi , he thinks that maybe a good solution could be to make use of std::functions when retrieving the material properties according to a certain TEMPERATURE instead of virtualizing the method. I am not really experienced in this issue so what do you guys think? @loumalouomega @philbucher

[loumalouomega]

After talking with @RiccardoRossi , he thinks that maybe a good solution could be to make use of std::functions when retrieving the material properties according to a certain TEMPERATURE instead of virtualizing the method. I am not really experienced in this issue so what do you guys think? @loumalouomega @philbucher

That's what I proposed...

[loumalouomega]

@AlejandroCornejo https://github.com/KratosMultiphysics/Kratos/pull/9053/commits/3028750d15a2378d7919dbd40c2154e90b269776 is an example, in this case using the function as a member variable (cons: takes more memory, one function for each GP). Only potential solution I can think, a template argument (variable prooperties for example), and then define a static function in compilation time, using std::conditional, otherwise I think is not possible to do it with a function

[RiccardoRossi]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call.

i'll expand soon about my idea on a possible use of the lambda functions

[loumalouomega]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call.

i'll expand soon about my idea on a possible use of the lambda functions

@RiccardoRossi you can see that I already implemented an initial version with lambda, the question is how to avoid the additional storage. The only way possible IMO is to add a template argument to the base laws and declare this lambda function as static

[loumalouomega]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call. i'll expand soon about my idea on a possible use of the lambda functions

@RiccardoRossi you can see that I already implemented an initial version with lambda, the question is how to avoid the additional storage. The only way possible IMO is to add a template argument to the base laws and declare this lambda function as static

See https://github.com/KratosMultiphysics/Kratos/commit/3028750d15a2378d7919dbd40c2154e90b269776

[philbucher]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call. i'll expand soon about my idea on a possible use of the lambda functions

@RiccardoRossi you can see that I already implemented an initial version with lambda, the question is how to avoid the additional storage. The only way possible IMO is to add a template argument to the base laws and declare this lambda function as static

See 3028750

not bad @loumalouomega

@RiccardoRossi we have to pay one price ... memory or runtime :) I would go for the solution of @loumalouomega with more memory, I don't like slowing down the computations

[loumalouomega]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call. i'll expand soon about my idea on a possible use of the lambda functions

@RiccardoRossi you can see that I already implemented an initial version with lambda, the question is how to avoid the additional storage. The only way possible IMO is to add a template argument to the base laws and declare this lambda function as static

See 3028750

not bad @loumalouomega

@RiccardoRossi we have to pay one price ... memory or runtime :) I would go for the solution of @loumalouomega with more memory, I don't like slowing down the computations

My second proposal does not slow down, nor consumes more memory (but will take longer to compile due to our beloved friend template arguments)

[philbucher]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call. i'll expand soon about my idea on a possible use of the lambda functions

@RiccardoRossi you can see that I already implemented an initial version with lambda, the question is how to avoid the additional storage. The only way possible IMO is to add a template argument to the base laws and declare this lambda function as static

See 3028750

not bad @loumalouomega @RiccardoRossi we have to pay one price ... memory or runtime :) I would go for the solution of @loumalouomega with more memory, I don't like slowing down the computations

My second proposal does not slow down, nor consumes more memory (but will take longer to compile due to our beloved friend template arguments)

but that one you haven't implemented right?

[loumalouomega]

ok, i just skimmed through the discussion and as a first comment please note that the benchmark should be done c++ only, otherwise what we will be measuring is the overhead of python insted of the overhead of the virtual call. i'll expand soon about my idea on a possible use of the lambda functions

@RiccardoRossi you can see that I already implemented an initial version with lambda, the question is how to avoid the additional storage. The only way possible IMO is to add a template argument to the base laws and declare this lambda function as static

See 3028750

not bad @loumalouomega @RiccardoRossi we have to pay one price ... memory or runtime :) I would go for the solution of @loumalouomega with more memory, I don't like slowing down the computations

My second proposal does not slow down, nor consumes more memory (but will take longer to compile due to our beloved friend template arguments)

but that one you haven't implemented right?

Nope, because @AlejandroCornejo has back-flashes from Vietnam every time a template argument is introduced in CL.

[philbucher]

Nope, because @AlejandroCornejo has back-flashes from Vietnam every time a template argument is introduced in CL.

he is not the only one!!!