[ShapeOpt] [WIP] Adding new response with stress aggregation for optimization

[dbaumgaertner]

This PR adds a new response function which allows for stress aggregation using the Kreisselmeier Steinhauser (KS) function.

By that, stress constraints may be considered over big parts of a structural domain whereas the given functionalities in Kratos only allowed for a single maximum value (whose location may heavily jump around). See figure below for an impression. There, a loaded hook structure is optimized for mass, while the stress in each element of the structure is not allowed to exceed a maximum value (dark red). The color plot shows the stresses on the surface for a constant value range. It can be seen, that the stresses significantly increase when reducing the mass (removing material). But they always stay within the specified limits.

This branch is not yet ready to be merged as it revealed an existing problem with thread safety in the current adjoint solution approach in case several elements are considered within the given stress response functions.

Also, this PR first requires an approval and merge of #4364.

Roadmap:

• Identify and fix the bug that destroys thread safety (remove race conditions)
• Remove temporary workaround with a manual set of the number of threads to 1
• Create a unit test for the stress aggregation response in the StructuralMechanicsApp
• Clean code and prepare for PR
• Wait for the merge of [StructuralMechanicsApplication][ShapeOpt] Added von mises stress sensitivities #4364

[dbaumgaertner]

Observations regarding the problematic race conditions so far:

1. A workaround that eliminates the race condition is to manually set the number of omp threads to 1 just before the StructuralMechanicsAdjointStaticSolver is initialized and reset it to the max values when it is finalized.

2. The problem gets visible in form of fluctuating shape sensitivities (dJdX) with reach new run.

3. The race condition only appears when a finite differencing is part of the computation of sensitivities. That is if all partial derivatives involving finite differencing (dpKdpX and dpJdpx) are hard set some values, then the shape sensitivity does not fluctuate anymore.

4. The adjoint variables are not subjected to the racing conditions and have constant values. This is because, in the case of stress constraints, the right-hand side of the adjoint system (-dpJdpu) is computed analytically and does not involve a finite differencing.

5. The fluctuating shape sensitivities appear because the function CalculateOnIntegrationPoints(VON_MISES_STRESS,...) sometimes returns meaningless or overly high stress values when it is called through the gradient functions within the class AdjointAggregatedStressResponseFunction. Noteworthy is that the stress values are fine when CalculateOnIntegrationPoints is called directly from a primal element (such as when the value of the response function is computed). By contrast, it returns nonsense values when called through the adjoint element which then internally calls the primal element (happens when the response function calculates its local gradients). So the problem seems to originate from the connection of primal and adjoint elements.

[philbucher]

minor

[philbucher]

blocking until this is resolved
As you pointed out this should not go into Master

[philbucher]

Suggested change

	elif response_type == "adjoint_nodal_displacement" or\
	elif response_type in ["adjoint_nodal_displacement", "adjoint_linear_strain_energy", ...]:

this would be another alternative of doing this

[armingeiser]

FYI @MFusseder the mentioned bug might be relevant for you as well.

[armingeiser]

The race condition only appears when a finite differencing is part of the computation of sensitivities. That is if all partial derivatives involving finite differencing (dpKdpX and dpJdpx) are hard set some values, then the shape sensitivity does not fluctuate anymore.

Is dpKdpX also causing problems? Because that is calculated in all other adjoint structural responses as well and there it seems to work (because it is called in serial).
dpJdpx until now was always just dependent on one node (or nodes of one element or even 0). Now with this response it is the first time that all nodal positions have an influence.

By contrast, it returns nonsense values when called through the adjoint element which then internally calls the primal element (happens when the response function calculates its local gradients). So the problem seems to originate from the connection of primal and adjoint elements.

The calculation of stress values at element level is threadsafe - so i would not expect that the calculation of simply the values does create any problems - also not when called from within the adjoints element at the primal element.

The finite differencing of derivatives w.r.t. node coordinates is not thread safe. Thats why i still guess that doing this in parallel is the root of the problem.

[MFusseder]

After using this new response for some days I have a comment. Within CalculateValue you compute also the values for the members mStressScalingFactor, mKSPrefactors and mSumPrefactors which are later needed in order to provide the partial derivatives. This approach perfectly works if you use structural_response.py. But if one uses directly a StructuralMechanicsAnalysis to solve the adjoint problem the function CalculateValueis not called. Therefor a extra custom running script would be needed. In principle I would propose to compute the values for member variables within InitializeSolutionStep of the response but this will lead to the problem that the primal model part is there not available. Right now I also have no better solution for this issue as to call CalculateValue but we should think about the fundamental problem that currently we have insufficient access to primal results (like stresses) during the computation of the partial derivatives which is needed for some responses as this example shows.

[dbaumgaertner]

As discussed in person: Indeed this is a problem. But for the sake of getting a first draft up and running, it was implemented this way. The main problem is, that the primal model part is only available in the CalculateValue() function. This again is due to the given structure of the adjoint analysis and requires a much deeper restructuring. Therefore, we first need to find a solution to the previous problem. That is why this PR is just WIP.