[ShapeOpt] [FluidDyn] Help needed solving a fluid ShapeOpt problem

[marcnunezc]

Hi! Im writing on behalf of the ExaQUte team! (tagging those who are in Kratos: @riccardotosi @RiccardoRossi @andreas-apostolatos and R. AMELA, Q. AYOUL–GUILMARD, S. GANESH, B. KEITH, A. KODAKKAL, M. RICCHIUTO )

This issue follows #4832, where I am using the drag sensitivities to solve a ShapeOpt problem for the same case. In this problem, I try to optimize a small box in 2D minimizing the drag force.

I have tried to run a small time analysis using steepest descent:

And penalized projection (with mass as a constraint):

I have some questions regarding the fluid dynamic case and the shape optimization app in general.

1. I see there is no fluid_response.py available (such as structural_response.py). Is there a file that I am missing out that I can use to solve fluid dynamics shape opt problems? For now I adapted the strucutral_response.py to the fluid case to do these first tests. Shall I open a PR with it? @sunethwarna @msandre

2. Regarding constraint functions, I can see that there is an option to set an "specified_value" and an equality/inequality condition. If I understand it correclty, the expected specified value is an specific value of the constraint function (for instance an specific value of mass = 1500). Could we add an option to set this specific value in terms of the initial one? For instance, specified_value = 0.8*initial_value

3. Also, related to the mass response function: I am using the mass as a constraint in the problem shown in #4832 (see .json below). The value of the mass seems to not change in my optimization analysis (I am putting a print here to check it), even though I have set an equality constraint with a differenet specified value than the initial. Is this expected? Unfortunately there are no examples in the test_examples folder in the ShapeOpt app that use the mass as a constraint that I can use to compare (unless I am missing out something). Using numbers, the mass value is always the initial: VALUE = 2484.000000000009 and the specifed value is: REFERENCE = 2500.0 I was expecting the mass value to get close to the specified value.

4. In the optimization output, I understand that the optimal design at each step is the one that results from the deformation of MESH_CHANGE (is what im using in the .gif above). Is this correct? Also, what is the difference between SHAPE_CHANGE and CONTROL_POINT_CHANGE?

5. Is there a way to get the primal analysis output (in this case VELOCITY/PRESSURE) of each optimizaiton step including mesh movement? I am using the gid_output_process in the FinalizeSolutionStep of my fluid_response but the mesh and the solution seem to be the same at every step. This makes me wonder if Im actually solving the deformed mesh or not at every step :disappointed:

6. Would you recommend an specific optimizer, filter settings, etc for this case? Providing that I am doing everything correctly, it seems that the optimizer does very few steps and the shape does not vary too much. This can also be related to point 3 where I am not seeing a mass change in the problem.

Finally I add here the parameters that I am using:

{
    "optimization_settings" : {
        "model_settings" : {
            "domain_size"           : 2,
            "model_part_name"       : "model",
            "model_import_settings" : {
                "input_type"     : "mdpa",
                "input_filename" : "ProblemZero"
            },
            "design_surface_sub_model_part_name" : "NoSlip2D_structure",
            "mesh_motion" : {
                "apply_mesh_solver" : true,
                "solver_settings" : {
                    "solver_type": "mesh_solver_structural_similarity",
                    "model_part_name"    : "model",
                    "model_import_settings"              : {
                        "input_type"     : "use_input_model_part"
                    },
                    "time_stepping" : {
                        "time_step"       : 1.0
                    },
                    "domain_size"     : 2,
                    "mesh_motion_linear_solver_settings"  : {
                        "solver_type"         : "amgcl",
                        "max_iteration"       : 500,
                        "tolerance"           : 1e-7,
                        "provide_coordinates" : false,
                        "smoother_type"       : "ilu0",
                        "krylov_type"         : "bicgstab",
                        "coarsening_type"     : "aggregation",
                        "scaling"             : false,
                        "verbosity"           : 1
                    },
                    "compute_reactions"     : false
                },
                "processes" : {
                    "boundary_conditions_process_list" : [{
                        "python_module" : "assign_vector_variable_process",
                        "kratos_module" : "KratosMultiphysics",
                        "help"          : "This process fixes the selected components of a given vector variable",
                        "process_name"  : "AssignVectorVariableProcess",
                        "Parameters"    : {
                            "mesh_id"         : 0,
                            "model_part_name" : "model.NoSlip2D_ground",
                            "variable_name"   : "MESH_DISPLACEMENT",
                            "value"           : [0.0,0.0,0.0]
                        }
                    },{
                        "python_module" : "assign_vector_variable_process",
                        "kratos_module" : "KratosMultiphysics",
                        "help"          : "This process fixes the selected components of a given vector variable",
                        "process_name"  : "AssignVectorVariableProcess",
                        "Parameters"    : {
                            "mesh_id"         : 0,
                            "model_part_name" : "model.AutomaticInlet2D_inlet",
                            "variable_name"   : "MESH_DISPLACEMENT",
                            "value"           : [0.0,0.0,0.0]
                        }
                    }]
                }
            }
        },
        "objectives" : [{
            "identifier" : "drag",
            "type"       : "minimization",
            "use_kratos" : true,
            "kratos_response_settings":{
                "primal_settings"  :"ProblemZeroParameters.json",
                "adjoint_settings" : "ProblemZeroAdjointParameters.json",
                "response_type"     : "drag",
                "model_part_name"       : "MainModelPart",
                "model_import_settings" : {
                    "input_type" : "use_input_model_part"
            },
            "gradient_mode"     : "semi_analytic",
            "step_size"         : 1e-9
            },
            "project_gradient_on_surface_normals" : true
        }],
        "constraints" : [{
            "identifier" : "mass",
            "type"       : "=",
            "reference"  : "specified_value",
            "reference_value" : 2500.0,
            "use_kratos" : true,
            "kratos_response_settings":{
                "response_type"          : "mass",
                "model_part_name" : "MainModelPart",
                "model_import_settings"              : {
                    "input_type"     : "use_input_model_part"
                },
                "material_import_settings" :{
                    "materials_filename": "materials_3D.json"
                },
                "primal_settings"        : "ProblemZeroParameters.json",
                "gradient_mode"          : "finite_differencing",
                "step_size"              : 1e-8
            },
            "project_gradient_on_surface_normals" : true
        }],
        "design_variables" : {
            "type"                : "vertex_morphing",
            "filter" : {
                "filter_function_type"        : "linear",
                "filter_radius"               : 2.0,
                "max_nodes_in_filter_radius"  : 10000,
                "matrix_free_filtering" : true
            }
        },
        "optimization_algorithm" : {
            "name"                    : "penalized_projection",
            "correction_scaling"      : 0.1,
            "use_adaptive_correction" : true,
            "max_iterations"          : 50,
            "relative_tolerance"      : 1e-3,
            "line_search" : {
                "line_search_type"           : "manual_stepping",
                "normalize_search_direction" : false,
                "step_size"                  : 1.0
            }
        },
        "output" : {
            "design_output_mode" : "WriteOptimizationModelPart",
            "nodal_results"      : [ "NORMALIZED_SURFACE_NORMAL",
                                     "DF1DX",
                                     "DC1DX",
                                     "DF1DX_MAPPED",
                                     "DC1DX_MAPPED",
                                     "CONTROL_POINT_UPDATE",
                                     "CONTROL_POINT_CHANGE",
                                     "SHAPE_UPDATE",
                                     "SHAPE_CHANGE",
                                     "MESH_CHANGE"],
            "output_format" : {
                "name": "gid"
            }
        }
    }
}

Thanks!!

[armingeiser]

I try to answer your points one by one:

1. There is no such file available yet. If you have a working state of it, you can of course make a MR, then we can discuss it directly on the code.

2. Thats currently not implemented (only maintaining the initial_value) - See e.g. the "test_examples/01_Hook..." A relative limit is of course something that could be implemented if needed.

3. This is not expected, at least you should see a slight change in the mass response value. Do i understand correctly that you are constraining the "mass" of the fluid domain?

4. SHAPE_CHANGE is the change of the mesh only at the design surface, MESH_CHANGE includes the mesh motion solution, where the SHAPE_CHANGE is propagated into the volume domain. CONTROL_POINT_CHANGE is the change of the control points of the optimization. They do NOT represent the actual geometry but are part of the Vertex Morphing theory (you can immagine them like control points of a NURBS geometry). You can neglect the CONTROL_POINT_CHANGE as it is an intermediate result.

5. If you use the same model_part for optimization and analysis, you should be able to switch on the MESH_CHANGE output in the analysis parameters. Then you can visualize the mesh change in GiD. Another option is to use individual model_parts for optimization and analysis. I think for more complex scenarios this is the better solution. I have done that for the case of the potential flow NACA optimization. Then you would get a individual folder for the primal/adjoint solution of each optimization step, also including the actual mesh aoutput to GiD. I can share the script if you like.

6. My suggestion is to do unconstrained optimization first until you understand what the optimization trys to do in this case. You can set the relative tolerance to a small or even negative value in order to force more optimization steps. I do not know the size f this cube, but you can also try to reduce the size of the filter radius. For the constraint example, you can try to increase the correction_scaling in order to get more contribution from the constraint when you are infeasible.

[marcnunezc]

@armingeiser

3.Yes, I wanted to let the square reduce its size a bit, therefore I configured the mass constraint with a specified value higher than the initial one.

5.Sure that would be great! Did you also create a potential_response for the potential case? I tried to do some tests with it in the branch cps/adding-potential-response.

6.Ok i will try your suggestions and get back to you!

[dbaumgaertner]

Just adding to the detailed explanations of @armingeiser:

3. From what I understood is, that you are actually measuring the mass of the fluid mesh. The way you specified the problem (target value higher than the initial value), the optimization should actually head towards the increased value. So there should definitely be a change in mass. Can you make sure, that the model part which the mass has, actually sees all these shape updates. A simple indirect test could be to run an unconstrained optimization (steepest descent) where you specify the mass as objective and let it minimize.

4. Generally, it is good advice, to first try an unconstrained case and then go for the constrained one. An option to consider in both cases is the line search. Here I would first set the option to normalize the search direction and then adjust the step size accordingly. So, specify normalize_search_direction = true and then think of the step size as the maximum nodal shape update per iteration in your geometric units. Meaning, if the cube has an edge length of 1, you could specify a step size of 0.25 to have reasonable updates in each step. Just note that this normalization, of course, might also lead to shape updates that are too big especially later in the optimization. To overcome this, you could change the "line_search_type" to "adaptive_stepping", which adjusts this value dynamically (the adaption only works for the unconstrained case). Anyways, these are parameters you could play with.

[marcnunezc]

Thanks @dbaumgaertner ,

As a quick update, after @armingeiser comments i realized that I was actually using two model_part_names ! (in the parameters I have "model" in some settings and "MainModelPart" in others), so that probably explains why I was seeing no change in the mass (since there were two model_parts).

However, right now when using one model_part for everything, the fluid solver (primal analysis) does not converge at the second step of the optimization. I am resetting the VELOCITY, ACCELERATION, PRESSURE, TIME and STEP of the primal analysis to the initial values after every optimization step. I guess there is still something i need to reset that I am not doing, but Im still looking into it!

[armingeiser]

@marcnunezc In this branch: feature/potential_flow_opt I have a test case for an optimization using the potential flow response applications/ShapeOptimizationApplication/test_examples/13_potential_flow_naca0012_real/run_optimization.py is the runscript that runs every potential flow simulation in a new folder.

This is a WIP, so of course there is room for improvement

[philbucher]

is this still open?

[marcnunezc]

We can close it