[Bug] Heat problems with unused conditions

[takeshimg92]

Description
TL;DR: keeping unused Conditions or SubModelParts on the MDPA for a heat problem yields wrong physics.
This is a problem I have noticed a few months back, only now reporting on it.

Considering the following simple heat diffusion problem in the ConvectionDiffusionApplication. I generate a tetrahedral mesh with 1263 nodes, corresponding to a beam which is 0.5m x 0.2m x 0.1cm along the x,y, z axes respectively.

Set up a Dirichlet condition at x=0 (surface x_1) for temperature T = 100 K. On the x=0.5 surface (call it x_2) set a flux of -25000 W/m2, corresponding to a total power output of 500W.

Setting thermal conductivity to 237 W/mK, the analytical solution for the temperature at x=0.5 is 47.257 K.

Scenario 1: construct MDPA with two submodelparts, x1 and x2, and only keep the Conditions that belong to these groups

• I get the correct value of 47.257 K to 3 decimal places.

Scenario 2: construct MDPA with more submodel parts (because why not?) -- create one at y=0 (y_1) and another at y=0.2 (y_2). Keep the Conditions that belong to these groups (but still only enforce boundary conditions on x_1 and x_2)

• Then, I get a wrong value of 37.930K at x=0.5.

Scenario 3: construct MDPA with just the two submodel parts x1 and x2, but keep all conditions (all the surfaces of the beam)

• Then, I get a wrong value of 16.3954 K at x=0.5

Scenario 4: construct MDPA with the 4 model parts above, but explicitly set the flux to zero at surfaces y_1 and y_2

• Then, I get a wrong value of 56.823 K at x=0.5.

I don't understand what is happening here.

Scope
ConvectionDiffusionApplication

To Reproduce

My solver settings are standard:

"solver_settings": {
                "model_part_name": "MyModelPart",
                "domain_size": 3,
                "echo_level": 1,
                "solver_type": 'stationary'
                "model_import_settings"   : {
			"input_type"                  : "mdpa",
			"input_filename"           : "mymesh"
                },
                "convection_diffusion_variables": {
                    "density_variable": "DENSITY",
                    "diffusion_variable": "CONDUCTIVITY",
                    "unknown_variable": "TEMPERATURE",
                    "volume_source_variable": "HEAT_FLUX",
                    "surface_source_variable": "FACE_HEAT_FLUX",
                    "projection_variable": "PROJECTED_SCALAR1",
                    "convection_variable": "CONVECTION_VELOCITY",
                    "mesh_velocity_variable": "MESH_VELOCITY",
                    "transfer_coefficient_variable": "",
                    "velocity_variable": "VELOCITY",
                    "specific_heat_variable": "SPECIFIC_HEAT",
                    "reaction_variable": "REACTION_FLUX",
                },
                "element_replace_settings": {
                    "element_name": "LaplacianElement",
                    "condition_name": "ThermalFace", 
                },
                "time_stepping": {"time_step": dt},
                "material_import_settings": {"materials_filename": model_part_name},
            },

I had the same error using ThermalFace and FluxCondition for element_replace_settings.

Expected behavior
In all cases, I would have expected temperature at x=0.5 to be equal to 47.257 K. Just having more conditions should not change the results if no additional boundary conditions are passed. Also, explicitly setting fluxes to zero at y_1, y_2 should be enough to replicate the original behavior.

Environment

• OS: Ubuntu 22.04
• Branch: master
• Python 3.10 and 3.12

[rubenzorrilla]

Hey @takeshimg92!

Thanks for the detailed report. At first glance, I think that the condition is adding some undesired contribution that results in something different to the do-nothing (adiabatic) that you would expect in this case.

This is related to the I/O and will most probably fixed once we fully adopt the geometry-based one, which will have full control on the entities that are substituted by ThermalFace thanks to the corresponding modeler.

If you still have them, it will be very handy for us to have the test cases corresponding to each of the scenarios.

Thanks in advance and sorry for the inconveniences.

[takeshimg92]

Hi @rubenzorrilla ! Thanks for the reply.
Can you help me better understand what you mean by this fully geometry based I/O? I think I am not up to date on this discussion.

[rubenzorrilla]

Sure. The point is that so far, we had two ways of providing the elements/conditions info:

1. To directly set them in the mdpa. This is for instance the approach followed in the StructuralMechanicsApplication, which normally require mixing several element types in one simulation.
2. To set "do-nothing" element and conditions (the base implementation indeed) in the mdpa and then substitute them alltogether by the corresponding type in the solver after reading the mdpa. This is what we did so far in applications such as the FluidDynamicsApplication or the ConvectionDiffusionApplication, which normally use a single element/condition type.

I suspect that the issue you're struggling with comes from the latter. Long story short, all conditions are being substituted by ThermalFace, including those that you just want to have in the mdpa for the sake of having them (something vary fair to be done).

The geometry-based I/O relies on only writing geometries in the mdpa, something that BTW is the unique thing we should ask to any preprocessor. In other workds, elements/conditions (i.e., physics) must be in the Kratos side. Then the geometries in each model part are substituted by the CreateEntitiesFromGeometries modeler according to the user needs. Note that these requires reading the mdpa at the AnalysisStage level. Then the input type in the solver_settings becomes the use_input_model_part, so we leave the responsability of setting up the mesh outside the solver.

[takeshimg92]

I see. Thank you for the explanation! For now, I will keep the issue open (if I have time I can push a branch with the test showing the problem).

[rubenzorrilla]

I see. Thank you for the explanation! For now, I will keep the issue open (if I have time I can push a branch with the test showing the problem).

I think that if you compressing and attachinc the case in here is enough (no need to create a test working with the testing framework).

[takeshimg92]

Hi @rubenzorrilla , I wrote a minimal example here: problematic_example.zip

As a note, in my original comment I was using a LaplacianElement; here I am using a MixedLaplacianElement because I wanted access to temperature gradients, but the conclusion is the same.

[rubenzorrilla]

Great. We'll try to have a look ASAP.