From 75dd85be20ad1f1e3834f4e392a2b3b73db7f241 Mon Sep 17 00:00:00 2001 From: maxnezddyur Date: Wed, 18 Dec 2024 13:37:57 -0700 Subject: [PATCH] add consistent scaling for optimization problems #29570 --- .../src/executioners/AdjointSolve.C | 29 ++- .../misc/scaling_test/forward_and_adjoint.i | 223 ++++++++++++++++++ .../misc/scaling_test/gold/main_out.json | 33 +++ .../test/tests/misc/scaling_test/main.i | 54 +++++ .../tests/misc/scaling_test/scaling_test.i | 123 ---------- .../test/tests/misc/scaling_test/tests | 34 ++- 6 files changed, 364 insertions(+), 132 deletions(-) create mode 100644 modules/optimization/test/tests/misc/scaling_test/forward_and_adjoint.i create mode 100644 modules/optimization/test/tests/misc/scaling_test/gold/main_out.json create mode 100644 modules/optimization/test/tests/misc/scaling_test/main.i delete mode 100644 modules/optimization/test/tests/misc/scaling_test/scaling_test.i diff --git a/modules/optimization/src/executioners/AdjointSolve.C b/modules/optimization/src/executioners/AdjointSolve.C index b49b55089223..053eaaae1eda 100644 --- a/modules/optimization/src/executioners/AdjointSolve.C +++ b/modules/optimization/src/executioners/AdjointSolve.C @@ -15,6 +15,7 @@ #include "NonlinearSystem.h" #include "NodalBCBase.h" +#include "libmesh/fuzzy_equals.h" #include "libmesh/petsc_matrix.h" #include "libmesh/petsc_vector.h" #include "petscmat.h" @@ -52,6 +53,14 @@ AdjointSolve::AdjointSolve(Executioner & ex) paramError("forward_system", "Forward system does not appear to be a 'NonlinearSystem'."); if (!dynamic_cast(&_nl_adjoint)) paramError("adjoint_system", "Adjoint system does not appear to be a 'NonlinearSystem'."); + // Adjoint system should never perform it's own automatic scaling. Scaling is + // taken from the forward system + _nl_adjoint.automaticScaling(false); + + // We need to force the forward system to have a scaling vector. This is + // incase a user provides scaling for indivial variables but doesn't have any + // AD objects. + _nl_forward.addScalingVector(); } bool @@ -95,6 +104,10 @@ AdjointSolve::solve() // Solve the adjoint system solver.adjoint_solve(matrix, solution, rhs, tol, maxits); + + // For scaling of the forward problem we need to apply correction factor + solution *= _nl_forward.getVector("scaling_factors"); + _nl_adjoint.update(); if (solver.get_converged_reason() < 0) { @@ -118,9 +131,6 @@ AdjointSolve::assembleAdjointSystem(SparseMatrix & matrix, const NumericVector & /*solution*/, NumericVector & rhs) { - if (_nl_adjoint.hasVector("scaling_factors")) - mooseWarning("Scaling factors are given to adjoint variables by the user. It is not necessary " - "to scale a adjoint system therefore the scaling factors will not be used."); _problem.computeJacobian(*_nl_forward.currentSolution(), matrix, _forward_sys_num); @@ -174,6 +184,19 @@ AdjointSolve::applyNodalBCs(SparseMatrix & matrix, void AdjointSolve::checkIntegrity() { + const auto adj_vars = _nl_adjoint.getVariables(0); + for (const auto & adj_var : adj_vars) + // If the user supplies any scaling factors for individual variables the + // adjoint system won't consistent. + if (!absolute_fuzzy_equals(adj_var->scalingFactor(), 1.0)) + mooseError("User supplied scaling factors for adjoint variables. Adjoint system is scaled " + "automatically by the forward system."); + + // This is to prevent automatic scaling of the adjoint system. Scaling is + // taken from the forward system + if (_nl_adjoint.hasVector("scaling_factors")) + _nl_adjoint.removeVector("scaling_factors"); + // Main thing is that the number of dofs in each system is the same if (_nl_forward.system().n_dofs() != _nl_adjoint.system().n_dofs()) mooseError( diff --git a/modules/optimization/test/tests/misc/scaling_test/forward_and_adjoint.i b/modules/optimization/test/tests/misc/scaling_test/forward_and_adjoint.i new file mode 100644 index 000000000000..eb15afb39e86 --- /dev/null +++ b/modules/optimization/test/tests/misc/scaling_test/forward_and_adjoint.i @@ -0,0 +1,223 @@ +[Mesh] + [generated] + type = GeneratedMeshGenerator + dim = 2 + nx = 10 + ny = 10 + ymin = 0 + xmin = 0 + xmax = 10 + ymax = 10 + [] +[] + +[Problem] + nl_sys_names = 'nl0 adjoint' + kernel_coverage_check = FALSE +[] + +[Variables] + [T_real] + initial_condition = 1e-8 + scaling = 10 + [] + [T_imag] + initial_condition = 1e-8 + [] + [T_real_adj] + solver_sys = adjoint + [] + [T_imag_adj] + solver_sys = adjoint + [] +[] + +[Kernels] + [heat_conduction_real] + type = MatDiffusion + variable = T_real + diffusivity = k + [] + [heat_source_real] + type = MatCoupledForce + variable = T_real + v = T_imag + material_properties = 'force_mat' + [] + + [heat_conduction_imag] + type = MatDiffusion + variable = T_imag + diffusivity = k + [] + [heat_source_imag] + type = MatCoupledForce + variable = T_imag + v = T_real + material_properties = 'force_mat' + coef = -1 + [] +[] + +[Materials] + [k_mat] + type = GenericFunctionMaterial + prop_names = 'k' + prop_values = 'kappa_func' + [] + [mats] + type = GenericConstantMaterial + prop_names = 'rho omega cp ' + prop_values = '1.0 1.0 1.0 ' + [] + [force_mat] + type = ParsedMaterial + property_name = force_mat + expression = 'rho * omega * cp' + material_property_names = 'rho omega cp' + [] + [phase] + type = ADParsedMaterial + coupled_variables = 'T_real T_imag' + expression = 'atan2(T_imag, T_real)' + property_name = phase + [] +[] + +[Functions] + [gauss] + type = ParsedFunction + expression = 'exp(-2.0 *(x^2 + y^2 + z^2)/(beam_radii^2))' + symbol_names = 'beam_radii' + symbol_values = '0.1' + [] + [kappa_func] + type = ParsedOptimizationFunction + expression = 'k ' + param_symbol_names = 'k ' + param_vector_name = 'params/k' + [] +[] + +[BCs] + [real_top] + type = FunctionNeumannBC + variable = T_real + boundary = top + function = 'exp((-2.0 *(x)^2)/0.1)' + [] +[] + +[DiracKernels] + [misfit_real] + type = ReporterPointSource + variable = T_real_adj + x_coord_name = measure_data_real/measurement_xcoord + y_coord_name = measure_data_real/measurement_ycoord + z_coord_name = measure_data_real/measurement_zcoord + value_name = measure_data_real/misfit_values + [] + [misfit_imag] + type = ReporterPointSource + variable = T_imag_adj + x_coord_name = measure_data_imag/measurement_xcoord + y_coord_name = measure_data_imag/measurement_ycoord + z_coord_name = measure_data_imag/measurement_zcoord + value_name = measure_data_imag/misfit_values + [] +[] + +[AuxVariables] + [phase] + [] +[] + +[AuxKernels] + [phase] + type = ParsedAux + variable = phase + coupled_variables = 'T_imag T_real' + expression = 'atan2(T_imag, T_real)' + execute_on = 'TIMESTEP_END' + [] +[] + +[Reporters] + [measure_data_real] + type = OptimizationData + variable = T_real + objective_name = objective_value + measurement_values = '0.10391 -0.0064' + measurement_points = '0.55 10 0 + 3.55 10 0' + [] + [measure_data_imag] + type = OptimizationData + objective_name = objective_value + variable = T_imag + measurement_values = '-0.08234 -0.00181' + measurement_points = '0.55 10 0 + 3.55 10 0' + [] + [params] + type = ConstantReporter + real_vector_names = 'k' + real_vector_values = '2' # Dummy value + [] + [gradient] + type = ParsedVectorReporter + name = inner + reporter_names = 'gradient_real/inner_product gradient_imag/inner_product' + reporter_symbols = 'a b' + expression = 'a+b' + execute_on = ADJOINT_TIMESTEP_END + execution_order_group = 1 + [] + [obj] + type = ParsedScalarReporter + name = value + reporter_names = 'measure_data_real/objective_value measure_data_imag/objective_value' + reporter_symbols = 'a b' + expression = 'a+b' + execute_on = ADJOINT_TIMESTEP_END + execution_order_group = 1 + [] +[] + +[VectorPostprocessors] + [gradient_real] + type = ElementOptimizationDiffusionCoefFunctionInnerProduct + variable = T_real_adj + forward_variable = T_real + function = kappa_func + execute_on = ADJOINT_TIMESTEP_END + [] + [gradient_imag] + type = ElementOptimizationDiffusionCoefFunctionInnerProduct + variable = T_imag_adj + forward_variable = T_imag + function = kappa_func + execute_on = ADJOINT_TIMESTEP_END + [] +[] + +[Executioner] + type = SteadyAndAdjoint + forward_system = nl0 + adjoint_system = adjoint + nl_rel_tol = 1e-12 + nl_abs_tol = 1e-10 + l_tol = 1e-10 + + automatic_scaling = false + + petsc_options_iname = '-pc_type' + petsc_options_value = 'lu' + line_search = l2 + verbose = true +[] + +[Outputs] + console = true + execute_on = FINAL +[] diff --git a/modules/optimization/test/tests/misc/scaling_test/gold/main_out.json b/modules/optimization/test/tests/misc/scaling_test/gold/main_out.json new file mode 100644 index 000000000000..86d94cb14786 --- /dev/null +++ b/modules/optimization/test/tests/misc/scaling_test/gold/main_out.json @@ -0,0 +1,33 @@ +{ + "reporters": { + "OptimizationReporter": { + "type": "GeneralOptimization", + "values": { + "grad_parameter_results": { + "type": "std::vector" + }, + "objective_value": { + "type": "double" + }, + "parameter_results": { + "type": "std::vector" + } + } + } + }, + "time_steps": [ + { + "OptimizationReporter": { + "grad_parameter_results": [ + -0.00141982365629666 + ], + "objective_value": 0.00025601754088785565, + "parameter_results": [ + 1.05 + ] + }, + "time": 0.0, + "time_step": 0 + } + ] +} diff --git a/modules/optimization/test/tests/misc/scaling_test/main.i b/modules/optimization/test/tests/misc/scaling_test/main.i new file mode 100644 index 000000000000..8e9f4a3ac2ef --- /dev/null +++ b/modules/optimization/test/tests/misc/scaling_test/main.i @@ -0,0 +1,54 @@ +[Optimization] +[] + +[OptimizationReporter] + type = GeneralOptimization + objective_name = objective_value + parameter_names = 'parameter_results' + num_values = '1' + initial_condition = '1.05' + lower_bounds = '0.001' +[] + +[Executioner] + type = Optimize + tao_solver = taobqnktr + petsc_options_iname = '-tao_max_it ' + petsc_options_value = '0 ' + verbose = true + output_optimization_iterations = true +[] + +[MultiApps] + [forward] + type = FullSolveMultiApp + input_files = forward_and_adjoint.i + execute_on = "FORWARD" + [] +[] + +[Transfers] + # FORWARD transfers + [toForward_measument] + type = MultiAppReporterTransfer + to_multi_app = forward + from_reporters = 'OptimizationReporter/parameter_results' + to_reporters = 'params/k' + [] + [fromForward] + type = MultiAppReporterTransfer + from_multi_app = forward + from_reporters = 'obj/value + gradient/inner' + to_reporters = 'OptimizationReporter/objective_value + OptimizationReporter/grad_parameter_results' + [] +[] + +[Outputs] + [out] + type = JSON + execute_on = 'FORWARD' + execute_system_information_on = NONE + [] +[] diff --git a/modules/optimization/test/tests/misc/scaling_test/scaling_test.i b/modules/optimization/test/tests/misc/scaling_test/scaling_test.i deleted file mode 100644 index 2394d166c0ed..000000000000 --- a/modules/optimization/test/tests/misc/scaling_test/scaling_test.i +++ /dev/null @@ -1,123 +0,0 @@ -[Problem] - nl_sys_names = 'nl0 adjoint' - kernel_coverage_check = false -[] - -[Mesh] - [gmg] - type = GeneratedMeshGenerator - dim = 3 - nx = 4 - ny = 2 - nz = 1 - xmin = 9.4615 - xmax = 92.0115 - ymin = 3.175 - ymax = 22.225 - zmin = 0.489 - zmax = 0.755 - [] -[] - -[Variables] - [T] - scaling = 10 - [] - [lam_T] - solver_sys = adjoint - scaling = 1e3 - [] -[] - -[Kernels] - [heat_conduction] - type = ADMatDiffusion - variable = T - diffusivity = thermal_conductivity - [] - [heat_source] - type = ADBodyForce - function = src_fuel_function - variable = T - [] -[] - -[BCs] - [dir_BC_front] - type = NeumannBC - variable = T - boundary = front - value = 2 - [] - [dir_BC_back] - type = DirichletBC - variable = T - boundary = back - value = 300 - [] -[] - -[Materials] - # fuel properties - [fuel_thermal] - type = ADGenericConstantMaterial - prop_names = thermal_conductivity - prop_values = 17.6e3 - [] -[] - -[Executioner] - type = SteadyAndAdjoint - forward_system = nl0 - adjoint_system = adjoint - - petsc_options_iname = '-pc_type -pc_factor_mat_solver_package' - petsc_options_value = 'lu mumps' - - line_search = 'none' - - nl_rel_tol = 1e-6 - nl_abs_tol = 1e-30 - nl_max_its = 10 - - l_max_its = 10 -[] - -##---------Forward Optimization stuff------------------# -[Reporters] - [measure_data] - type = ConstantReporter - real_vector_names = 'x y z u weight' - real_vector_values = '0.2 0.2 0.0; 0.3 0.8 0.0; 0 0 0; 5 5 5; 1 1 1' - [] - [params_fuel] - type = ConstantReporter - real_vector_names = 'source' - real_vector_values = '5e7' # Dummy - [] -[] - -[Functions] - [src_fuel_function] - type = ParsedOptimizationFunction - expression = q - param_symbol_names = 'q' - param_vector_name = 'params_fuel/source' - [] -[] -##---------Adjoint Optimization stuff------------------# -[DiracKernels] - [adjointLoad_T] - type = ReporterPointSource - variable = lam_T - x_coord_name = measure_data/x - y_coord_name = measure_data/y - z_coord_name = measure_data/z - value_name = measure_data/u - [] -[] -##--------- Outputs ------------------# - -[Debug] - show_var_residual_norms = true -[] diff --git a/modules/optimization/test/tests/misc/scaling_test/tests b/modules/optimization/test/tests/misc/scaling_test/tests index 4ac2ba534afd..b4aae457ec74 100644 --- a/modules/optimization/test/tests/misc/scaling_test/tests +++ b/modules/optimization/test/tests/misc/scaling_test/tests @@ -1,10 +1,32 @@ [Tests] - issues = '#25646' + issues = '#25646 #29570' design = 'SteadyAndAdjoint.md' - [scaling_factor_for_adjoint_variables] - type = 'RunException' - input = scaling_test.i - expect_err = 'Scaling factors are given to adjoint variables by the user. It is not necessary to scale a adjoint system therefore the scaling factors will not be used.' - requirement = "The system shall throw out a warning when scalings are assigned to adjoint variable to remind users that no scaling will be performed for adjoint system." + [optimization_scaling] + requirement = "The system shall be able to scale forward and adjoint systems consistently by" + + [manual_scaling] + type = JSONDiff + input = main.i + jsondiff = main_out.json + detail = "by manually providing scaling factors for the forward problem variables," + [] + [automatic_scaling] + type = JSONDiff + input = main.i + jsondiff = main_out.json + cli_args = "MultiApps/forward/cli_args='Executioner/automatic_scaling=true'" + detail = "or by using the automatic scaling option." + [] + [] + [optimization_scaling_errors] + [scaling_factor_for_adjoint_variables] + type = 'RunException' + input = main.i + cli_args = "MultiApps/forward/cli_args='Variables/T_real_adj/scaling=5'" + expect_err = 'User supplied scaling factors for adjoint variables. Adjoint system is scaled automatically by the forward system.' + requirement = "The system shall throw out an error when scalings are manually assigned to " + "adjoint variables to remind users that no scaling will be performed for adjoint " + "system." + [] [] []