diff --git a/applications/ChimeraApplication/python_scripts/navier_stokes_solver_fractionalstep_chimera.py b/applications/ChimeraApplication/python_scripts/navier_stokes_solver_fractionalstep_chimera.py index 5800cbf6994a..7bff884ef94a 100755 --- a/applications/ChimeraApplication/python_scripts/navier_stokes_solver_fractionalstep_chimera.py +++ b/applications/ChimeraApplication/python_scripts/navier_stokes_solver_fractionalstep_chimera.py @@ -19,7 +19,7 @@ class NavierStokesSolverFractionalStepForChimera(NavierStokesSolverFractionalSte def __init__(self, model, custom_settings): [self.chimera_settings, self.chimera_internal_parts, custom_settings] = chimera_setup_utils.SeparateAndValidateChimeraSettings(custom_settings) super(NavierStokesSolverFractionalStepForChimera,self).__init__(model,custom_settings) - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStepForChimera", "Construction of NavierStokesSolverFractionalStepForChimera finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesSolverFractionalStepForChimera finished.") def AddVariables(self): super(NavierStokesSolverFractionalStepForChimera,self).AddVariables() @@ -31,9 +31,7 @@ def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosChimera.ROTATION_MESH_DISPLACEMENT) self.main_model_part.AddNodalSolutionStepVariable(KratosChimera.ROTATION_MESH_VELOCITY) - - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStepForChimera", "Fluid solver variables added correctly.") - + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid chimera solver variables added correctly.") def ImportModelPart(self): if(self.settings["model_import_settings"]["input_type"].GetString() == "chimera"): @@ -45,63 +43,25 @@ def ImportModelPart(self): material_file_name = self.settings["material_import_settings"]["materials_filename"].GetString() import KratosMultiphysics.ChimeraApplication.chimera_modelpart_import as chim_mp_imp chim_mp_imp.ImportChimeraModelparts(self.main_model_part, chimera_mp_import_settings, material_file=material_file_name, parallel_type="OpenMP") - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStepForChimera", " Import of all chimera modelparts completed.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, " Import of all chimera modelparts completed.") else:# we can use the default implementation in the base class super(NavierStokesSolverFractionalStepForChimera,self).ImportModelPart() def Initialize(self): - self.chimera_process = chimera_setup_utils.GetApplyChimeraProcess(self.model, self.chimera_settings, self.settings) - self.computing_model_part =self.main_model_part - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._get_automatic_time_stepping_utility() - - #TODO: next part would be much cleaner if we passed directly the parameters to the c++ - if self.settings["consider_periodic_conditions"] == True: - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStepForChimera Periodic conditions are not implemented in this case .") - raise NotImplementedError - else: - self.solver_settings = KratosChimera.FractionalStepSettings(self.computing_model_part, - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.settings["time_order"].GetInt(), - self.settings["use_slip_conditions"].GetBool(), - self.settings["move_mesh_flag"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool()) - - self.solver_settings.SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.solver_settings.SetStrategy(KratosCFD.StrategyLabel.Velocity, - self.velocity_linear_solver, - self.settings["velocity_tolerance"].GetDouble(), - self.settings["maximum_velocity_iterations"].GetInt()) - - self.solver_settings.SetStrategy(KratosCFD.StrategyLabel.Pressure, - self.pressure_linear_solver, - self.settings["pressure_tolerance"].GetDouble(), - self.settings["maximum_pressure_iterations"].GetInt()) - - - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.solver = KratosCFD.FSStrategy(self.computing_model_part, - self.solver_settings, - self.settings["predictor_corrector"].GetBool(), - KratosCFD.PATCH_INDEX) - else: - self.solver = KratosChimera.FSStrategyForChimera(self.computing_model_part, - self.solver_settings, - self.settings["predictor_corrector"].GetBool()) - - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) + # Call the base solver to create the solution strategy + super(NavierStokesSolverFractionalStepForChimera,self).Initialize() - (self.solver).Initialize() - - self.solver.Check() - - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStepForChimera", "Solver initialization finished.") - - chimera_setup_utils.SetChimeraInternalPartsFlag(self.model, self.chimera_internal_parts) + # Chimera utilities initialization + self.chimera_process = chimera_setup_utils.GetApplyChimeraProcess( + self.model, + self.chimera_settings, + self.settings) + chimera_setup_utils.SetChimeraInternalPartsFlag( + self.model, + self.chimera_internal_parts) + def GetComputingModelPart(self): + return self.main_model_part def InitializeSolutionStep(self): self.chimera_process.ExecuteInitializeSolutionStep() @@ -110,4 +70,53 @@ def InitializeSolutionStep(self): def FinalizeSolutionStep(self): super(NavierStokesSolverFractionalStepForChimera,self).FinalizeSolutionStep() ## Depending on the setting this will clear the created constraints - self.chimera_process.ExecuteFinalizeSolutionStep() \ No newline at end of file + self.chimera_process.ExecuteFinalizeSolutionStep() + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + domain_size = computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + + # Create the pressure and velocity linear solvers + # Note that linear_solvers is a tuple. The first item is the pressure + # linear solver. The second item is the velocity linear solver. + linear_solvers = self._GetLinearSolver() + + # Create the fractional step settings instance + # TODO: next part would be much cleaner if we passed directly the parameters to the c++ + if self.settings["consider_periodic_conditions"].GetBool(): + KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStepForChimera Periodic conditions are not implemented in this case .") + raise NotImplementedError + else: + fractional_step_settings = KratosChimera.FractionalStepSettings( + computing_model_part, + domain_size, + self.settings["time_order"].GetInt(), + self.settings["use_slip_conditions"].GetBool(), + self.settings["move_mesh_flag"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool()) + + # Set the strategy echo level + fractional_step_settings.SetEchoLevel(self.settings["echo_level"].GetInt()) + + # Set the velocity and pressure fractional step strategy settings + fractional_step_settings.SetStrategy(KratosCFD.StrategyLabel.Pressure, + linear_solvers[0], + self.settings["pressure_tolerance"].GetDouble(), + self.settings["maximum_pressure_iterations"].GetInt()) + + fractional_step_settings.SetStrategy(KratosCFD.StrategyLabel.Velocity, + linear_solvers[1], + self.settings["velocity_tolerance"].GetDouble(), + self.settings["maximum_velocity_iterations"].GetInt()) + + # Create the fractional step strategy + if self.settings["consider_periodic_conditions"].GetBool() == True: + KratosMultiphysics.Logger.PrintInfo("FSStrategyForChimera Periodic conditions are not implemented in this case .") + raise NotImplementedError + else: + solution_strategy = KratosChimera.FSStrategyForChimera( + computing_model_part, + fractional_step_settings, + self.settings["predictor_corrector"].GetBool()) + + return solution_strategy \ No newline at end of file diff --git a/applications/ChimeraApplication/python_scripts/navier_stokes_solver_vmsmonolithic_chimera.py b/applications/ChimeraApplication/python_scripts/navier_stokes_solver_vmsmonolithic_chimera.py index 008dd6f96f56..60ae029067f9 100644 --- a/applications/ChimeraApplication/python_scripts/navier_stokes_solver_vmsmonolithic_chimera.py +++ b/applications/ChimeraApplication/python_scripts/navier_stokes_solver_vmsmonolithic_chimera.py @@ -4,8 +4,6 @@ import KratosMultiphysics import KratosMultiphysics.ChimeraApplication as KratosChimera from KratosMultiphysics.ChimeraApplication import chimera_setup_utils -# Import applications -import KratosMultiphysics.FluidDynamicsApplication as KratosCFD # Import base class file from KratosMultiphysics.FluidDynamicsApplication.navier_stokes_solver_vmsmonolithic import NavierStokesSolverMonolithic @@ -17,7 +15,7 @@ class NavierStokesSolverMonolithicChimera(NavierStokesSolverMonolithic): def __init__(self, model, custom_settings): [self.chimera_settings, self.chimera_internal_parts, custom_settings] = chimera_setup_utils.SeparateAndValidateChimeraSettings(custom_settings) super(NavierStokesSolverMonolithicChimera,self).__init__(model,custom_settings) - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicChimera", "Construction of NavierStokesSolverMonolithic finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesSolverMonolithicChimera finished.") def AddVariables(self): super(NavierStokesSolverMonolithicChimera,self).AddVariables() @@ -28,7 +26,7 @@ def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosChimera.ROTATION_MESH_DISPLACEMENT) self.main_model_part.AddNodalSolutionStepVariable(KratosChimera.ROTATION_MESH_VELOCITY) - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicChimera", "Fluid solver variables added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid chimera solver variables added correctly.") def ImportModelPart(self): if(self.settings["model_import_settings"]["input_type"].GetString() == "chimera"): @@ -40,106 +38,40 @@ def ImportModelPart(self): material_file_name = self.settings["material_import_settings"]["materials_filename"].GetString() import KratosMultiphysics.ChimeraApplication.chimera_modelpart_import as chim_mp_imp chim_mp_imp.ImportChimeraModelparts(self.main_model_part, chimera_mp_import_settings, material_file=material_file_name, parallel_type="OpenMP") - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicChimera", " Import of all chimera modelparts completed.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, " Import of all chimera modelparts completed.") else:# we can use the default implementation in the base class super(NavierStokesSolverMonolithicChimera,self).ImportModelPart() def Initialize(self): - self.chimera_process = chimera_setup_utils.GetApplyChimeraProcess(self.model, self.chimera_settings, self.settings) - self.computing_model_part = self.main_model_part - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._get_automatic_time_stepping_utility() + # Call the base solver to create the solution strategy + super(NavierStokesSolverMonolithicChimera,self).Initialize() - # Creating the solution strategy - self.conv_criteria = KratosCFD.VelPrCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) - - # Creating the time integration scheme - if (self.element_integrates_in_time): - # "Fake" scheme for those cases in where the element manages the time integration - # It is required to perform the nodal update once the current time step is solved - self.time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip( - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) - # In case the BDF2 scheme is used inside the element, the BDF time discretization utility is required to update the BDF coefficients - if (self.settings["time_scheme"].GetString() == "bdf2"): - time_order = 2 - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) - else: - err_msg = "Requested elemental time scheme \"" + self.settings["time_scheme"].GetString()+ "\" is not available.\n" - err_msg += "Available options are: \"bdf2\"" - raise Exception(err_msg) - else: - if not hasattr(self, "_turbulence_model_solver"): - # Bossak time integration scheme - if self.settings["time_scheme"].GetString() == "bossak": - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.time_scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( - self.settings["alpha"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - KratosCFD.PATCH_INDEX) - else: - self.time_scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( - self.settings["alpha"].GetDouble(), - self.settings["move_mesh_strategy"].GetInt(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - # BDF2 time integration scheme - elif self.settings["time_scheme"].GetString() == "bdf2": - self.time_scheme = KratosCFD.GearScheme() - # Time scheme for steady state fluid solver - elif self.settings["time_scheme"].GetString() == "steady": - self.time_scheme = KratosCFD.ResidualBasedSimpleSteadyScheme( - self.settings["velocity_relaxation"].GetDouble(), - self.settings["pressure_relaxation"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - else: - err_msg = "Requested time scheme " + self.settings["time_scheme"].GetString() + " is not available.\n" - err_msg += "Available options are: \"bossak\", \"bdf2\" and \"steady\"" - raise Exception(err_msg) - else: - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicForChimera turbulent solver is not possible.") - raise NotImplementedError - - if self.settings["consider_periodic_conditions"].GetBool() == True: - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicForChimera Periodic conditions are not implemented in this case .") - raise NotImplementedError - else: - builder_and_solver = KratosChimera.ResidualBasedBlockBuilderAndSolverWithConstraintsForChimera(self.linear_solver) - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(self.computing_model_part, - self.time_scheme, - self.linear_solver, - self.conv_criteria, - builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.formulation.SetProcessInfo(self.computing_model_part) - - (self.solver).Initialize() - - self.solver.Check() - - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicChimera", "Solver initialization finished.") - - chimera_setup_utils.SetChimeraInternalPartsFlag(self.model, self.chimera_internal_parts) + # Chimera utilities initialization + self.chimera_process = chimera_setup_utils.GetApplyChimeraProcess( + self.model, + self.chimera_settings, + self.settings) + chimera_setup_utils.SetChimeraInternalPartsFlag( + self.model, + self.chimera_internal_parts) + def GetComputingModelPart(self): + return self.main_model_part def InitializeSolutionStep(self): self.chimera_process.ExecuteInitializeSolutionStep() super(NavierStokesSolverMonolithicChimera,self).InitializeSolutionStep() - def FinalizeSolutionStep(self): super(NavierStokesSolverMonolithicChimera,self).FinalizeSolutionStep() ## Depending on the setting this will clear the created constraints self.chimera_process.ExecuteFinalizeSolutionStep() + + def _CreateBuilderAndSolver(self): + linear_solver = self._GetLinearSolver() + if self.settings["consider_periodic_conditions"].GetBool(): + KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithicForChimera Periodic conditions are not implemented in this case .") + raise NotImplementedError + else: + builder_and_solver = KratosChimera.ResidualBasedBlockBuilderAndSolverWithConstraintsForChimera(linear_solver) + return builder_and_solver diff --git a/applications/ChimeraApplication/tests/chimera_analysis_base_test.py b/applications/ChimeraApplication/tests/chimera_analysis_base_test.py index ca14b006e189..be035d1ced1e 100644 --- a/applications/ChimeraApplication/tests/chimera_analysis_base_test.py +++ b/applications/ChimeraApplication/tests/chimera_analysis_base_test.py @@ -22,7 +22,29 @@ def _run_test(self,settings_file_name): "process_name" : "VtkOutputProcess", "help" : "This process writes postprocessing files for Paraview", "Parameters" : { - "model_part_name" : "FluidModelPart", + "model_part_name" : "FluidModelPart.Parts_background_surface", + "output_control_type" : "step", + "output_frequency" : 1, + "file_format" : "ascii", + "output_precision" : 3, + "output_sub_model_parts" : false, + "write_deformed_configuration" : true, + "folder_name" : "test_vtk_output", + "save_output_files_in_folder" : true, + "nodal_solution_step_data_variables" : ["VELOCITY","PRESSURE","DISTANCE","MESH_VELOCITY"], + "nodal_data_value_variables" : [], + "element_flags" : ["ACTIVE"], + "nodal_flags" : ["VISITED","CHIMERA_INTERNAL_BOUNDARY"], + "element_data_value_variables" : [], + "condition_data_value_variables" : [] + } + },{ + "python_module" : "vtk_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "VtkOutputProcess", + "help" : "This process writes postprocessing files for Paraview", + "Parameters" : { + "model_part_name" : "FluidModelPart.Parts_patch_surface", "output_control_type" : "step", "output_frequency" : 1, "file_format" : "ascii", diff --git a/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_adjoint_solver.py b/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_adjoint_solver.py index 7c8b4d0bd011..9accf30a2034 100644 --- a/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_adjoint_solver.py +++ b/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_adjoint_solver.py @@ -75,7 +75,7 @@ def __init__(self, model, custom_settings): self.element_name = self.formulation.element_name self.condition_name = self.formulation.condition_name - KratosMultiphysics.Logger.PrintInfo("::[PotentialFlowAdjointSolver]:: ", "Construction finished") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction finished") def AddVariables(self): super(PotentialFlowAdjointSolver, self).AddVariables() @@ -84,50 +84,67 @@ def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.SHAPE_SENSITIVITY) self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.NORMAL_SENSITIVITY) - KratosMultiphysics.Logger.PrintInfo("::[PotentialFlowAdjointSolver]:: ", "Variables ADDED") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Variables ADDED") def AddDofs(self): KratosMultiphysics.VariableUtils().AddDof(KCPFApp.ADJOINT_VELOCITY_POTENTIAL, self.main_model_part) KratosMultiphysics.VariableUtils().AddDof(KCPFApp.ADJOINT_AUXILIARY_VELOCITY_POTENTIAL, self.main_model_part) def Initialize(self): - self._ComputeNodalNeighbours() + # Call base solver Initialize() to calculate the nodal neighbours and initialize the strategy + super(PotentialFlowAdjointSolver, self).Initialize() - """Perform initialization after adding nodal variables and dofs to the main model part. """ - if self.response_function_settings["response_type"].GetString() == "adjoint_lift_jump_coordinates": - self.response_function = KCPFApp.AdjointLiftJumpCoordinatesResponseFunction(self.main_model_part, self.response_function_settings) - else: - raise Exception("invalid response_type: " + self.response_function_settings["response_type"].GetString()) - - self.sensitivity_builder=KratosMultiphysics.SensitivityBuilder(self.sensitivity_settings,self.main_model_part, self.response_function) - self.sensitivity_builder.Initialize() - - scheme = KratosMultiphysics.ResidualBasedAdjointStaticScheme(self.response_function) - - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - self.solver = KratosMultiphysics.ResidualBasedLinearStrategy( - self.main_model_part, - scheme, - self.linear_solver, - builder_and_solver, - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["calculate_solution_norm"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) + # Initialize the response function and the sensitivity builder + self._GetResponseFunction().Initialize() + self._GetSensitivityBuilder().Initialize() - self.solver.SetEchoLevel(self.settings["echo_level"].GetInt()) - self.solver.Check() - - self.response_function.Initialize() - - KratosMultiphysics.Logger.PrintInfo("::[PotentialFlowAdjointSolver]:: ", "Finished initialization.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Finished initialization.") def InitializeSolutionStep(self): super(PotentialFlowAdjointSolver, self).InitializeSolutionStep() - self.response_function.InitializeSolutionStep() + self._GetResponseFunction().InitializeSolutionStep() def FinalizeSolutionStep(self): super(PotentialFlowAdjointSolver, self).FinalizeSolutionStep() - self.response_function.FinalizeSolutionStep() - self.sensitivity_builder.UpdateSensitivities() - + self._GetResponseFunction().FinalizeSolutionStep() + self._GetSensitivityBuilder().UpdateSensitivities() + + @classmethod + def _GetStrategyType(self): + strategy_type = "linear" + return strategy_type + + def _CreateScheme(self): + # Fake scheme creation to do the solution update + response_function = self._GetResponseFunction() + scheme = KratosMultiphysics.ResidualBasedAdjointStaticScheme(response_function) + return scheme + + def _GetResponseFunction(self): + if not hasattr(self, '_response_function'): + self._response_function = self.__CreateResponseFunction() + return self._response_function + + def __CreateResponseFunction(self): + computing_model_part = self.GetComputingModelPart() + if self.response_function_settings["response_type"].GetString() == "adjoint_lift_jump_coordinates": + response_function = KCPFApp.AdjointLiftJumpCoordinatesResponseFunction( + computing_model_part, + self.response_function_settings) + else: + raise Exception("Invalid response_type: " + self.response_function_settings["response_type"].GetString()) + return response_function + + def _GetSensitivityBuilder(self): + if not hasattr(self, '_sensitivity_builder'): + self._sensitivity_builder = self.__CreateSensitivityBuilder() + return self._sensitivity_builder + + def __CreateSensitivityBuilder(self): + computing_model_part = self.GetComputingModelPart() + response_function = self._GetResponseFunction() + sensitivity_builder = KratosMultiphysics.SensitivityBuilder( + self.sensitivity_settings, + computing_model_part, + response_function) + return sensitivity_builder diff --git a/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_solver.py b/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_solver.py index 45c5fc029fa2..06fe9646ab37 100644 --- a/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_solver.py +++ b/applications/CompressiblePotentialFlowApplication/python_scripts/potential_flow_solver.py @@ -1,5 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.CompressiblePotentialFlowApplication as KCPFApp @@ -140,9 +138,6 @@ def __init__(self, model, custom_settings): self._validate_settings_in_baseclass=True # To be removed eventually super(PotentialFlowSolver, self).__init__(model, custom_settings) - # There is only a single rank in OpenMP, we always print - self._is_printing_rank = True - # Set the element and condition names for the replace settings self.formulation = PotentialFlowFormulation(self.settings["formulation"]) self.element_name = self.formulation.element_name @@ -154,10 +149,6 @@ def __init__(self, model, custom_settings): self.main_model_part.ProcessInfo.SetValue(KCPFApp.REFERENCE_CHORD,self.reference_chord) self.element_has_nodal_properties = False - #construct the linear solvers - import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory - self.linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - def AddVariables(self): # Degrees of freedom self.main_model_part.AddNodalSolutionStepVariable(KCPFApp.VELOCITY_POTENTIAL) @@ -169,67 +160,84 @@ def AddVariables(self): variable = KratosMultiphysics.KratosGlobals.GetVariable(variable_name) self.main_model_part.AddNodalSolutionStepVariable(variable) - KratosMultiphysics.Logger.PrintInfo("::[PotentialFlowSolver]:: ", "Variables ADDED") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Variables ADDED") def AddDofs(self): KratosMultiphysics.VariableUtils().AddDof(KCPFApp.VELOCITY_POTENTIAL, self.main_model_part) KratosMultiphysics.VariableUtils().AddDof(KCPFApp.AUXILIARY_VELOCITY_POTENTIAL, self.main_model_part) def Initialize(self): - self._ComputeNodalNeighbours() - - time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticScheme() - strategy = self._GetStrategyType() - if strategy == "linear": - # TODO: Rename to self.strategy once we upgrade the base FluidDynamicsApplication solvers - self.solver = KratosMultiphysics.ResidualBasedLinearStrategy( - self.GetComputingModelPart(), - time_scheme, - self.linear_solver, - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["calculate_solution_norm"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - elif strategy == "non_linear": - conv_criteria = KratosMultiphysics.ResidualCriteria( - self.settings["relative_tolerance"].GetDouble(), - self.settings["absolute_tolerance"].GetDouble()) - max_iterations = self.settings["maximum_iterations"].GetInt() - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy( - self.GetComputingModelPart(), - time_scheme, - self.linear_solver, - conv_criteria, - max_iterations, - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - else: - raise Exception("Element not implemented") + self._ComputeNodalElementalNeighbours() - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - self.solver.Initialize() + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) + solution_strategy.Initialize() + + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def AdvanceInTime(self, current_time): raise Exception("AdvanceInTime is not implemented. Potential Flow simulations are steady state.") - def _ComputeNodalNeighbours(self): + def _ComputeNodalElementalNeighbours(self): # Find nodal neigbours util call - KratosMultiphysics.FindNodalNeighboursProcess(self.main_model_part).Execute() + data_communicator = KratosMultiphysics.DataCommunicator.GetDefault() + KratosMultiphysics.FindGlobalNodalElementalNeighboursProcess( + data_communicator, + self.GetComputingModelPart()).Execute() def _GetStrategyType(self): element_type = self.settings["formulation"]["element_type"].GetString() if "incompressible" in element_type: if not self.settings["formulation"].Has("penalty_coefficient"): - strategy = "linear" + strategy_type = "linear" elif self.settings["formulation"]["penalty_coefficient"].GetDouble() == 0.0: - strategy = "linear" + strategy_type = "linear" else: - strategy = "non_linear" + strategy_type = "non_linear" elif "compressible" in element_type: - strategy = "non_linear" + strategy_type = "non_linear" else: - strategy = "" + strategy_type = None + return strategy_type - return strategy + @classmethod + def _CreateScheme(self): + # Fake scheme creation to do the solution update + scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticScheme() + return scheme + + def _CreateConvergenceCriterion(self): + convergence_criterion = KratosMultiphysics.ResidualCriteria( + self.settings["relative_tolerance"].GetDouble(), + self.settings["absolute_tolerance"].GetDouble()) + return convergence_criterion + + def _CreateSolutionStrategy(self): + strategy_type = self._GetStrategyType() + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + if strategy_type == "linear": + solution_strategy = KratosMultiphysics.ResidualBasedLinearStrategy( + computing_model_part, + time_scheme, + linear_solver, + self.settings["compute_reactions"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + self.settings["calculate_solution_norm"].GetBool(), + self.settings["move_mesh_flag"].GetBool()) + elif strategy_type == "non_linear": + convergence_criterion = self._GetConvergenceCriterion() + solution_strategy = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy( + computing_model_part, + time_scheme, + linear_solver, + convergence_criterion, + self.settings["maximum_iterations"].GetInt(), + self.settings["compute_reactions"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + self.settings["move_mesh_flag"].GetBool()) + else: + err_msg = "Unknown strategy type: \'" + strategy_type + "\'. Valid options are \'linear\' and \'non_linear\'." + raise Exception(err_msg) + return solution_strategy \ No newline at end of file diff --git a/applications/CompressiblePotentialFlowApplication/tests/embedded_test/embedded_circle_adjoint_parameters.json b/applications/CompressiblePotentialFlowApplication/tests/embedded_test/embedded_circle_adjoint_parameters.json index 0ecb0217e041..f38bc38b6ea5 100644 --- a/applications/CompressiblePotentialFlowApplication/tests/embedded_test/embedded_circle_adjoint_parameters.json +++ b/applications/CompressiblePotentialFlowApplication/tests/embedded_test/embedded_circle_adjoint_parameters.json @@ -23,7 +23,6 @@ "step_size" : 1e-6 }, "sensitivity_settings" : { - "sensitivity_model_part_name" : "Parts_Parts_Auto1", "nodal_solution_step_sensitivity_variables" : ["NORMAL_SENSITIVITY"], "nodal_solution_step_sensitivity_calculation_is_thread_safe" : false, "build_mode": "static" diff --git a/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_analytical_parameters.json b/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_analytical_parameters.json index 4c3cd96b5a6f..d7003fd8a522 100644 --- a/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_analytical_parameters.json +++ b/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_analytical_parameters.json @@ -23,7 +23,6 @@ "gradient_mode" : "analytic" }, "sensitivity_settings" : { - "sensitivity_model_part_name" : "Parts_Parts_Auto1", "nodal_solution_step_sensitivity_variables" : ["SHAPE_SENSITIVITY"], "nodal_solution_step_sensitivity_calculation_is_thread_safe" : false, "build_mode": "static" diff --git a/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_parameters.json b/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_parameters.json index ead526436170..c02f7d65b4ce 100644 --- a/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_parameters.json +++ b/applications/CompressiblePotentialFlowApplication/tests/naca0012_small_adjoint_test/naca0012_small_sensitivities_adjoint_parameters.json @@ -24,7 +24,6 @@ "step_size" : 1e-9 }, "sensitivity_settings" : { - "sensitivity_model_part_name" : "Parts_Parts_Auto1", "nodal_solution_step_sensitivity_variables" : ["SHAPE_SENSITIVITY"], "nodal_solution_step_sensitivity_calculation_is_thread_safe" : false, "build_mode": "static" @@ -83,4 +82,4 @@ }], "auxiliar_process_list" :[] } -} \ No newline at end of file +} diff --git a/applications/FSIApplication/python_scripts/partitioned_embedded_fsi_base_solver.py b/applications/FSIApplication/python_scripts/partitioned_embedded_fsi_base_solver.py index 50d767f6a9a1..3da8a0505620 100644 --- a/applications/FSIApplication/python_scripts/partitioned_embedded_fsi_base_solver.py +++ b/applications/FSIApplication/python_scripts/partitioned_embedded_fsi_base_solver.py @@ -172,13 +172,13 @@ def Predict(self): self.__UpdateLevelSet() # Correct the updated level set - self.fluid_solver._GetDistanceModificationProcess().ExecuteInitializeSolutionStep() + self.fluid_solver.GetDistanceModificationProcess().ExecuteInitializeSolutionStep() # Fluid solver prediction self.fluid_solver.Predict() # Restore the fluid node fixity to its original status - self.fluid_solver._GetDistanceModificationProcess().ExecuteFinalizeSolutionStep() + self.fluid_solver.GetDistanceModificationProcess().ExecuteFinalizeSolutionStep() def GetComputingModelPart(self): err_msg = 'Calling GetComputingModelPart() method in a partitioned solver.\n' diff --git a/applications/FluidDynamicsApplication/python_scripts/adjoint_fluid_solver.py b/applications/FluidDynamicsApplication/python_scripts/adjoint_fluid_solver.py index f4040e416c96..5efdc4baa7be 100755 --- a/applications/FluidDynamicsApplication/python_scripts/adjoint_fluid_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/adjoint_fluid_solver.py @@ -1,45 +1,22 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 -import sys - # Importing the Kratos Library import KratosMultiphysics -from KratosMultiphysics.python_solver import PythonSolver +# from KratosMultiphysics.python_solver import PythonSolver import KratosMultiphysics.FluidDynamicsApplication as KratosCFD -from KratosMultiphysics.FluidDynamicsApplication import check_and_prepare_model_process_fluid +from KratosMultiphysics.FluidDynamicsApplication.fluid_solver import FluidSolver def CreateSolver(model, custom_settings): return AdjointFluidSolver(model, custom_settings) -class AdjointFluidSolver(PythonSolver): +class AdjointFluidSolver(FluidSolver): def __init__(self, model, settings): - super(AdjointFluidSolver,self).__init__(model, settings) - ## Set the element and condition names for the replace settings - ## These should be defined in derived classes - self.element_name = None - self.condition_name = None + # Overwrite the default buffer size in base FluidSolver + # TODO: CHECK WHY THE DEFAULT BUFFER SIZE IS 2 IF THE DEFAULT TIME SCHEME IS BOSSAK self.min_buffer_size = 2 - # Either retrieve the model part from the model or create a new one - model_part_name = self.settings["model_part_name"].GetString() - - if model_part_name == "": - raise Exception('Please specify a model_part name!') - - if self.model.HasModelPart(model_part_name): - self.main_model_part = self.model.GetModelPart(model_part_name) - else: - self.main_model_part = self.model.CreateModelPart(model_part_name) - - domain_size = self.settings["domain_size"].GetInt() - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DOMAIN_SIZE, domain_size) - - def AddVariables(self): - raise Exception("Trying to call AdjointFluidSolver.AddVariables(). Implement the AddVariables() method in the specific derived solver.") - def AddDofs(self): KratosMultiphysics.VariableUtils().AddDof(KratosCFD.ADJOINT_FLUID_VECTOR_1_X, self.main_model_part) KratosMultiphysics.VariableUtils().AddDof(KratosCFD.ADJOINT_FLUID_VECTOR_1_Y, self.main_model_part) @@ -48,80 +25,33 @@ def AddDofs(self): KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Adjoint fluid solver DOFs added correctly.") - def ImportModelPart(self): - # we can use the default implementation in the base class - self._ImportModelPart(self.main_model_part,self.settings["model_import_settings"]) - - def PrepareModelPart(self): - if not self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED]: - ## Set fluid properties from materials json file - materials_imported = self._SetPhysicalProperties() - if not materials_imported: - KratosMultiphysics.Logger.PrintWarning(self.__class__.__name__, "Material properties have not been imported. Check \'material_import_settings\' in your ProjectParameters.json.") - ## Replace default elements and conditions - self._ReplaceElementsAndConditions() - ## Executes the check and prepare model process - self._ExecuteCheckAndPrepare() - ## Set buffer size - self.main_model_part.SetBufferSize(self.min_buffer_size) - - KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Model reading finished.") - - def ExportModelPart(self): - ## Model part writing - name_out_file = self.settings["model_import_settings"]["input_filename"].GetString()+".out" - KratosMultiphysics.ModelPartIO(name_out_file, KratosMultiphysics.IO.WRITE).WriteModelPart(self.main_model_part) - - KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Model export finished.") - - def GetMinimumBufferSize(self): - return self.min_buffer_size - - def Initialize(self): - raise Exception("Calling AdjointFluidSolver.Initialize() base method. Please implement a custom Initialize() method for your solver.") - - def AdvanceInTime(self, current_time): - dt = self._ComputeDeltaTime() - new_time = current_time + dt - - self.main_model_part.CloneTimeStep(new_time) - self.main_model_part.ProcessInfo[KratosMultiphysics.STEP] += 1 - - return new_time - def InitializeSolutionStep(self): - self.solver.InitializeSolutionStep() - self.response_function.InitializeSolutionStep() + self._GetSolutionStrategy().InitializeSolutionStep() + self.GetResponseFunction().InitializeSolutionStep() if hasattr(self, "_adjoint_turbulence_model_solver"): self._adjoint_turbulence_model_solver.InitializeSolutionStep() def Predict(self): - self.solver.Predict() + self._GetSolutionStrategy().Predict() def SolveSolutionStep(self): - return self.solver.SolveSolutionStep() + return self._GetSolutionStrategy().SolveSolutionStep() def FinalizeSolutionStep(self): - (self.solver).FinalizeSolutionStep() - self.response_function.FinalizeSolutionStep() + self._GetSolutionStrategy().FinalizeSolutionStep() + self.GetResponseFunction().FinalizeSolutionStep() if hasattr(self, "_adjoint_turbulence_model_solver"): self._adjoint_turbulence_model_solver.FinalizeSolutionStep() - self.sensitivity_builder.UpdateSensitivities() + self.GetSensitivityBuilder().UpdateSensitivities() def Check(self): - (self.solver).Check() + self._GetSolutionStrategy().Check() if hasattr(self, "_adjoint_turbulence_model_solver"): self._adjoint_turbulence_model_solver.Check() - def Clear(self): - (self.solver).Clear() - - def GetComputingModelPart(self): - return self.main_model_part.GetSubModelPart("fluid_computational_model_part") - def _ReplaceElementsAndConditions(self): ## Get number of nodes and domain size elem_num_nodes = self._GetElementNumNodes() @@ -136,6 +66,7 @@ def _ReplaceElementsAndConditions(self): cond_num_nodes = domain_size ## Complete the element name + # TODO: EXPORT THE ADJOINT FOLLOWING THE CONVENTION. ONCE THIS IS DONE WE CAN USE THE FUNCTION IN THE BASE CLASS if (self.element_name is not None): new_elem_name = self.element_name + str(int(domain_size)) + "D" else: @@ -155,42 +86,6 @@ def _ReplaceElementsAndConditions(self): ## Call the replace elements and conditions process KratosMultiphysics.ReplaceElementsAndConditionsProcess(self.main_model_part, self.settings["element_replace_settings"]).Execute() - def _GetElementNumNodes(self): - if self.main_model_part.NumberOfElements() != 0: - if sys.version_info[0] >= 3: # python3 syntax - element_num_nodes = len(self.main_model_part.Elements.__iter__().__next__().GetNodes()) - else: # python2 syntax - element_num_nodes = len(self.main_model_part.Elements.__iter__().next().GetNodes()) - else: - element_num_nodes = 0 - - element_num_nodes = self.main_model_part.GetCommunicator().GetDataCommunicator().MaxAll(element_num_nodes) - return element_num_nodes - - def _GetConditionNumNodes(self): - if self.main_model_part.NumberOfConditions() != 0: - if sys.version_info[0] >= 3: # python3 syntax - condition_num_nodes = len(self.main_model_part.Conditions.__iter__().__next__().GetNodes()) - else: # python2 syntax - condition_num_nodes = len(self.main_model_part.Conditions.__iter__().next().GetNodes()) - else: - condition_num_nodes = 0 - - condition_num_nodes = self.main_model_part.GetCommunicator().GetDataCommunicator().MaxAll(condition_num_nodes) - return condition_num_nodes - - def _ExecuteCheckAndPrepare(self): - ## Check that the input read has the shape we like - prepare_model_part_settings = KratosMultiphysics.Parameters("{}") - prepare_model_part_settings.AddValue("volume_model_part_name",self.settings["volume_model_part_name"]) - prepare_model_part_settings.AddValue("skin_parts",self.settings["skin_parts"]) - - check_and_prepare_model_process_fluid.CheckAndPrepareModelProcess(self.main_model_part, prepare_model_part_settings).Execute() - - current_buffer_size = self.main_model_part.GetBufferSize() - if(self.GetMinimumBufferSize() > current_buffer_size): - self.main_model_part.SetBufferSize( self.GetMinimumBufferSize() ) - def _ComputeDeltaTime(self): if self.settings["time_stepping"]["automatic_time_step"].GetBool(): raise Exception("Automatic time stepping is not supported by adjoint fluid solver.") @@ -198,24 +93,6 @@ def _ComputeDeltaTime(self): delta_time = self.settings["time_stepping"]["time_step"].GetDouble() return delta_time - def _SetPhysicalProperties(self): - # Check if the fluid properties are provided using a .json file - materials_filename = self.settings["material_import_settings"]["materials_filename"].GetString() - if (materials_filename != ""): - # Add constitutive laws and material properties from json file to model parts. - material_settings = KratosMultiphysics.Parameters("""{"Parameters": {"materials_filename": ""}} """) - material_settings["Parameters"]["materials_filename"].SetString(materials_filename) - KratosMultiphysics.ReadMaterialsUtility(material_settings, self.model) - materials_imported = True - else: - materials_imported = False - - # If the element uses nodal material properties, transfer them to the nodes - if self.element_has_nodal_properties: - self._SetNodalProperties() - - return materials_imported - def _SetNodalProperties(self): # Get density and dynamic viscostity from the properties of the first element for el in self.main_model_part.Elements: @@ -230,6 +107,45 @@ def _SetNodalProperties(self): else: raise Exception("No fluid elements found in the main model part.") # Transfer the obtained properties to the nodes - KratosMultiphysics.VariableUtils().SetScalarVar(KratosMultiphysics.DENSITY, rho, self.main_model_part.Nodes) + KratosMultiphysics.VariableUtils().SetVariable(KratosMultiphysics.DENSITY, rho, self.main_model_part.Nodes) + # In RANS full adjoints, viscosity is summation of kinematic viscosity and turubulent viscosity. + # turbulent viscosity is read from hdf5, and viscosity also must be read from hdf5. + # therefore viscosity should not be updated only with kinematic viscosity. if not hasattr(self, "_adjoint_turbulence_model_solver"): - KratosMultiphysics.VariableUtils().SetScalarVar(KratosMultiphysics.VISCOSITY, kin_viscosity, self.main_model_part.Nodes) + KratosMultiphysics.VariableUtils().SetVariable(KratosMultiphysics.VISCOSITY, kin_viscosity, self.main_model_part.Nodes) + + def GetResponseFunction(self): + if not hasattr(self, '_response_function'): + self._response_function = self.__CreateResponseFunction() + return self._response_function + + def __CreateResponseFunction(self): + domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + response_type = self.settings["response_function_settings"]["response_type"].GetString() + if response_type == "drag": + if domain_size == 2: + response_function = KratosCFD.DragResponseFunction2D( + self.settings["response_function_settings"]["custom_settings"], + self.main_model_part) + elif domain_size == 3: + response_function = KratosCFD.DragResponseFunction3D( + self.settings["response_function_settings"]["custom_settings"], + self.main_model_part) + else: + raise Exception("Invalid DOMAIN_SIZE: " + str(domain_size)) + else: + raise Exception("Invalid response_type: " + response_type + ". Available response functions: \'drag\'.") + return response_function + + def GetSensitivityBuilder(self): + if not hasattr(self, '_sensitivity_builder'): + self._sensitivity_builder = self.__CreateSensitivityBuilder() + return self._sensitivity_builder + + def __CreateSensitivityBuilder(self): + response_function = self.GetResponseFunction() + sensitivity_builder = KratosMultiphysics.SensitivityBuilder( + self.settings["sensitivity_settings"], + self.main_model_part, + response_function) + return sensitivity_builder diff --git a/applications/FluidDynamicsApplication/python_scripts/adjoint_vmsmonolithic_solver.py b/applications/FluidDynamicsApplication/python_scripts/adjoint_vmsmonolithic_solver.py index 9954556f599b..4a972dd1c0c1 100644 --- a/applications/FluidDynamicsApplication/python_scripts/adjoint_vmsmonolithic_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/adjoint_vmsmonolithic_solver.py @@ -1,4 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 # Importing the Kratos Library import KratosMultiphysics @@ -51,7 +50,8 @@ def GetDefaultSettings(cls): "time_step" : -0.1 }, "adjoint_turbulence_model_solver_settings": {}, - "consider_periodic_conditions": false + "consider_periodic_conditions": false, + "assign_neighbour_elements_to_conditions": false }""") default_settings.AddMissingParameters(super(AdjointVMSMonolithicSolver, cls).GetDefaultSettings()) @@ -62,10 +62,6 @@ def __init__(self, model, custom_settings): super(AdjointVMSMonolithicSolver,self).__init__(model,custom_settings) self.element_has_nodal_properties = True - # construct the linear solver - import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory - self.linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - if not self.settings["adjoint_turbulence_model_solver_settings"].IsEquivalentTo(KratosMultiphysics.Parameters("{}")): # if not empty self._adjoint_turbulence_model_solver = CreateAdjointTurbulenceModel(model, self.settings["adjoint_turbulence_model_solver_settings"]) @@ -78,11 +74,12 @@ def __init__(self, model, custom_settings): elif self.settings["domain_size"].GetInt() == 3: self.condition_name = "SurfaceCondition" + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of AdjointVMSMonolithicSolver finished.") def AddVariables(self): - self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.VELOCITY) self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.ACCELERATION) self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.PRESSURE) @@ -99,6 +96,7 @@ def AddVariables(self): # Adding variables required for the turbulence modelling if hasattr(self, "_adjoint_turbulence_model_solver"): + # TODO: THIS HAS TO BE A METHOD SetFluidModelPart() IN THE ADJOINT TURBULENCE MODEL SOLVER self._adjoint_turbulence_model_solver.fluid_model_part = self.main_model_part self._adjoint_turbulence_model_solver.AddVariables() @@ -114,58 +112,49 @@ def AddDofs(self): self._adjoint_turbulence_model_solver.AddDofs() def Initialize(self): + # Construct and set the solution strategy + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) - self.computing_model_part = self.GetComputingModelPart() - - self.response_function = self.GetResponseFunction() - - self.sensitivity_builder = KratosMultiphysics.SensitivityBuilder(self.settings["sensitivity_settings"], self.main_model_part, self.response_function) - - if self.settings["scheme_settings"]["scheme_type"].GetString() == "bossak": - self.time_scheme = KratosMultiphysics.ResidualBasedAdjointBossakScheme(self.settings["scheme_settings"], self.response_function) - elif self.settings["scheme_settings"]["scheme_type"].GetString() == "steady": - self.time_scheme = KratosMultiphysics.ResidualBasedAdjointSteadyScheme(self.response_function) - else: - raise Exception("invalid scheme_type: " + self.settings["scheme_settings"]["scheme_type"].GetString()) - - if self.settings["consider_periodic_conditions"].GetBool() == True: - builder_and_solver = KratosCFD.ResidualBasedBlockBuilderAndSolverPeriodic( - self.linear_solver, KratosCFD.PATCH_INDEX) - else: - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - - self.solver = KratosMultiphysics.ResidualBasedLinearStrategy(self.main_model_part, - self.time_scheme, - self.linear_solver, - builder_and_solver, - False, - False, - False, - False) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) - + # If there is adjoint turbulence model, initialize it if hasattr(self, "_adjoint_turbulence_model_solver"): self._adjoint_turbulence_model_solver.Initialize() - (self.solver).Initialize() - (self.response_function).Initialize() - (self.sensitivity_builder).Initialize() - KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") + # Initialize the strategy and adjoint utilities + solution_strategy.Initialize() + self.GetResponseFunction().Initialize() + self.GetSensitivityBuilder().Initialize() - def GetResponseFunction(self): - domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") - if self.settings["response_function_settings"]["response_type"].GetString() == "drag": - if (domain_size == 2): - return KratosCFD.DragResponseFunction2D(self.settings["response_function_settings"]["custom_settings"], self.main_model_part) - elif (domain_size == 3): - return KratosCFD.DragResponseFunction3D(self.settings["response_function_settings"]["custom_settings"], self.main_model_part) - else: - raise Exception("Invalid DOMAIN_SIZE: " + str(domain_size)) + def _CreateScheme(self): + response_function = self.GetResponseFunction() + scheme_type = self.settings["scheme_settings"]["scheme_type"].GetString() + if scheme_type == "bossak": + scheme = KratosMultiphysics.ResidualBasedAdjointBossakScheme( + self.settings["scheme_settings"], + response_function) + elif scheme_type == "steady": + scheme = KratosMultiphysics.ResidualBasedAdjointSteadyScheme(response_function) else: - raise Exception("invalid response_type: " + self.settings["response_function_settings"]["response_type"].GetString()) \ No newline at end of file + raise Exception("Invalid scheme_type: " + scheme_type) + return scheme + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + builder_and_solver = self._GetBuilderAndSolver() + calculate_reaction_flag = False + reform_dof_set_at_each_step = False + calculate_norm_dx_flag = False + move_mesh_flag = False + return KratosMultiphysics.ResidualBasedLinearStrategy( + computing_model_part, + time_scheme, + linear_solver, + builder_and_solver, + calculate_reaction_flag, + reform_dof_set_at_each_step, + calculate_norm_dx_flag, + move_mesh_flag) diff --git a/applications/FluidDynamicsApplication/python_scripts/fluid_solver.py b/applications/FluidDynamicsApplication/python_scripts/fluid_solver.py index b5a60fdd0f17..08bc20bcc768 100755 --- a/applications/FluidDynamicsApplication/python_scripts/fluid_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/fluid_solver.py @@ -1,10 +1,9 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - import sys # Importing the Kratos Library import KratosMultiphysics from KratosMultiphysics.python_solver import PythonSolver +import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory # Import applications import KratosMultiphysics.FluidDynamicsApplication as KratosCFD @@ -14,7 +13,30 @@ def CreateSolver(model, custom_settings): return FluidSolver(model, custom_settings) class FluidSolver(PythonSolver): + """The base class for fluid dynamics solvers. + + This class provides functions for importing and exporting models, + adding nodal variables and dofs and solving each solution step. + + Depending on the formulation type, derived classes may require to + override some (or all) the following functions: + + _CreateScheme + _CreateConvergenceCriterion + _CreateLinearSolver + _CreateBuilderAndSolver + _CreateSolutionStrategy + The solution strategy, builder_and_solver, etc. should alway be retrieved + using the getter functions _GetSolutionStrategy, _GetBuilderAndSolver, + etc. from this base class. + + Only the member variables listed below should be accessed directly. + + Public member variables: + model -- the model containing the modelpart used to construct the solver. + settings -- Kratos parameters containing solver settings. + """ def __init__(self, model, settings): super(FluidSolver,self).__init__(model, settings) @@ -51,14 +73,14 @@ def AddDofs(self): KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.VELOCITY_Z, KratosMultiphysics.REACTION_Z,self.main_model_part) KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.PRESSURE, KratosMultiphysics.REACTION_WATER_PRESSURE,self.main_model_part) - KratosMultiphysics.Logger.PrintInfo("FluidSolver", "Fluid solver DOFs added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid solver DOFs added correctly.") def ImportModelPart(self): # we can use the default implementation in the base class self._ImportModelPart(self.main_model_part,self.settings["model_import_settings"]) def PrepareModelPart(self): - if not self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED]: + if not self.is_restarted(): ## Set fluid properties from materials json file materials_imported = self._SetPhysicalProperties() if not materials_imported: @@ -70,14 +92,14 @@ def PrepareModelPart(self): ## Set buffer size self.main_model_part.SetBufferSize(self.min_buffer_size) - KratosMultiphysics.Logger.PrintInfo("FluidSolver", "Model reading finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Model reading finished.") def ExportModelPart(self): ## Model part writing name_out_file = self.settings["model_import_settings"]["input_filename"].GetString()+".out" KratosMultiphysics.ModelPartIO(name_out_file, KratosMultiphysics.IO.WRITE).WriteModelPart(self.main_model_part) - KratosMultiphysics.Logger.PrintInfo("FluidSolver", "Model export finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Model export finished.") def GetMinimumBufferSize(self): return self.min_buffer_size @@ -96,32 +118,32 @@ def AdvanceInTime(self, current_time): def InitializeSolutionStep(self): if self._TimeBufferIsInitialized(): - self.solver.InitializeSolutionStep() + self._GetSolutionStrategy().InitializeSolutionStep() def Predict(self): if self._TimeBufferIsInitialized(): - self.solver.Predict() + self._GetSolutionStrategy().Predict() def SolveSolutionStep(self): if self._TimeBufferIsInitialized(): - is_converged = self.solver.SolveSolutionStep() + is_converged = self._GetSolutionStrategy().SolveSolutionStep() if not is_converged: msg = "Fluid solver did not converge for step " + str(self.main_model_part.ProcessInfo[KratosMultiphysics.STEP]) + "\n" msg += "corresponding to time " + str(self.main_model_part.ProcessInfo[KratosMultiphysics.TIME]) + "\n" - KratosMultiphysics.Logger.PrintWarning("FluidSolver",msg) + KratosMultiphysics.Logger.PrintWarning(self.__class__.__name__, msg) return is_converged else: return True def FinalizeSolutionStep(self): if self._TimeBufferIsInitialized(): - (self.solver).FinalizeSolutionStep() + self._GetSolutionStrategy().FinalizeSolutionStep() def Check(self): - (self.solver).Check() + self._GetSolutionStrategy().Check() def Clear(self): - (self.solver).Clear() + self._GetSolutionStrategy().Clear() def GetComputingModelPart(self): if not self.main_model_part.HasSubModelPart("fluid_computational_model_part"): @@ -208,23 +230,13 @@ def _ExecuteCheckAndPrepare(self): def _ComputeDeltaTime(self): # Automatic time step computation according to user defined CFL number if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - delta_time = self.EstimateDeltaTimeUtility.EstimateDt() + delta_time = self.GetEstimateDtUtility().EstimateDt() # User-defined delta time else: delta_time = self.settings["time_stepping"]["time_step"].GetDouble() return delta_time - def _GetAutomaticTimeSteppingUtility(self): - if (self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2): - EstimateDeltaTimeUtility = KratosCFD.EstimateDtUtility2D(self.GetComputingModelPart(), - self.settings["time_stepping"]) - else: - EstimateDeltaTimeUtility = KratosCFD.EstimateDtUtility3D(self.GetComputingModelPart(), - self.settings["time_stepping"]) - - return EstimateDeltaTimeUtility - def _SetPhysicalProperties(self): # Check if the fluid properties are provided using a .json file materials_filename = self.settings["material_import_settings"]["materials_filename"].GetString() @@ -259,3 +271,162 @@ def _SetNodalProperties(self): # Transfer the obtained properties to the nodes KratosMultiphysics.VariableUtils().SetVariable(KratosMultiphysics.DENSITY, rho, self.main_model_part.Nodes) KratosMultiphysics.VariableUtils().SetVariable(KratosMultiphysics.VISCOSITY, kin_viscosity, self.main_model_part.Nodes) + + # TODO: I THINK THIS SHOULD BE MOVED TO THE BASE PYTHON SOLVER + def is_restarted(self): + # this function avoids the long call to ProcessInfo and is also safer + # in case the detection of a restart is changed later + return self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED] + + def GetEstimateDtUtility(self): + if not hasattr(self, '_estimate_dt_utility'): + self._estimate_dt_utility = self._CreateEstimateDtUtility() + return self._estimate_dt_utility + + def _GetScheme(self): + if not hasattr(self, '_scheme'): + self._scheme = self._CreateScheme() + return self._scheme + + def _GetConvergenceCriterion(self): + if not hasattr(self, '_convergence_criterion'): + self._convergence_criterion = self._CreateConvergenceCriterion() + return self._convergence_criterion + + def _GetLinearSolver(self): + if not hasattr(self, '_linear_solver'): + self._linear_solver = self._CreateLinearSolver() + return self._linear_solver + + def _GetBuilderAndSolver(self): + if not hasattr(self, '_builder_and_solver'): + self._builder_and_solver = self._CreateBuilderAndSolver() + return self._builder_and_solver + + def _GetSolutionStrategy(self): + if not hasattr(self, '_solution_strategy'): + self._solution_strategy = self._CreateSolutionStrategy() + return self._solution_strategy + + def _CreateEstimateDtUtility(self): + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + if domain_size == 2: + estimate_dt_utility = KratosCFD.EstimateDtUtility2D( + self.GetComputingModelPart(), + self.settings["time_stepping"]) + else: + estimate_dt_utility = KratosCFD.EstimateDtUtility3D( + self.GetComputingModelPart(), + self.settings["time_stepping"]) + + return estimate_dt_utility + + def _CreateScheme(self): + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + # Cases in which the element manages the time integration + if self.element_integrates_in_time: + # "Fake" scheme for those cases in where the element manages the time integration + # It is required to perform the nodal update once the current time step is solved + scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip( + domain_size, + domain_size + 1) + # In case the BDF2 scheme is used inside the element, the BDF time discretization utility is required to update the BDF coefficients + if (self.settings["time_scheme"].GetString() == "bdf2"): + time_order = 2 + self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) + else: + err_msg = "Requested elemental time scheme \"" + self.settings["time_scheme"].GetString()+ "\" is not available.\n" + err_msg += "Available options are: \"bdf2\"" + raise Exception(err_msg) + # Cases in which a time scheme manages the time integration + else: + # Time scheme without turbulence modelling + if not hasattr(self, "_turbulence_model_solver"): + # Bossak time integration scheme + if self.settings["time_scheme"].GetString() == "bossak": + if self.settings["consider_periodic_conditions"].GetBool() == True: + scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( + self.settings["alpha"].GetDouble(), + domain_size, + KratosCFD.PATCH_INDEX) + else: + scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( + self.settings["alpha"].GetDouble(), + self.settings["move_mesh_strategy"].GetInt(), + domain_size) + # BDF2 time integration scheme + elif self.settings["time_scheme"].GetString() == "bdf2": + scheme = KratosCFD.GearScheme() + # Time scheme for steady state fluid solver + elif self.settings["time_scheme"].GetString() == "steady": + scheme = KratosCFD.ResidualBasedSimpleSteadyScheme( + self.settings["velocity_relaxation"].GetDouble(), + self.settings["pressure_relaxation"].GetDouble(), + domain_size) + else: + err_msg = "Requested time scheme " + self.settings["time_scheme"].GetString() + " is not available.\n" + err_msg += "Available options are: \"bossak\", \"bdf2\" and \"steady\"" + raise Exception(err_msg) + # Time scheme with turbulence modelling + else: + self._turbulence_model_solver.Initialize() + if self.settings["time_scheme"].GetString() == "bossak": + scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( + self.settings["alpha"].GetDouble(), + self.settings["move_mesh_strategy"].GetInt(), + domain_size, + self.settings["turbulence_model_solver_settings"]["velocity_pressure_relaxation_factor"].GetDouble(), + self._turbulence_model_solver.GetTurbulenceSolvingProcess()) + # Time scheme for steady state fluid solver + elif self.settings["time_scheme"].GetString() == "steady": + scheme = KratosCFD.ResidualBasedSimpleSteadyScheme( + self.settings["velocity_relaxation"].GetDouble(), + self.settings["pressure_relaxation"].GetDouble(), + domain_size, + self._turbulence_model_solver.GetTurbulenceSolvingProcess()) + return scheme + + def _CreateLinearSolver(self): + linear_solver_configuration = self.settings["linear_solver_settings"] + return linear_solver_factory.ConstructSolver(linear_solver_configuration) + + def _CreateConvergenceCriterion(self): + if self.settings["time_scheme"].GetString() == "steady": + convergence_criterion = KratosMultiphysics.ResidualCriteria( + self.settings["relative_velocity_tolerance"].GetDouble(), + self.settings["absolute_velocity_tolerance"].GetDouble()) + else: + convergence_criterion = KratosCFD.VelPrCriteria( + self.settings["relative_velocity_tolerance"].GetDouble(), + self.settings["absolute_velocity_tolerance"].GetDouble(), + self.settings["relative_pressure_tolerance"].GetDouble(), + self.settings["absolute_pressure_tolerance"].GetDouble()) + convergence_criterion.SetEchoLevel(self.settings["echo_level"].GetInt()) + return convergence_criterion + + def _CreateBuilderAndSolver(self): + linear_solver = self._GetLinearSolver() + if self.settings["consider_periodic_conditions"].GetBool(): + builder_and_solver = KratosCFD.ResidualBasedBlockBuilderAndSolverPeriodic( + linear_solver, + KratosCFD.PATCH_INDEX) + else: + builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(linear_solver) + return builder_and_solver + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + convergence_criterion = self._GetConvergenceCriterion() + builder_and_solver = self._GetBuilderAndSolver() + return KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy( + computing_model_part, + time_scheme, + linear_solver, + convergence_criterion, + builder_and_solver, + self.settings["maximum_iterations"].GetInt(), + self.settings["compute_reactions"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + self.settings["move_mesh_flag"].GetBool()) diff --git a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_ale_fluid_solver.py b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_ale_fluid_solver.py index ef37d0f26592..a37b6a3c8ccf 100644 --- a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_ale_fluid_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_ale_fluid_solver.py @@ -1,5 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics as KM @@ -29,11 +27,11 @@ def _ManipulateFluidSolverSettingsForReactionsComputation(cls, fluid_solver_sett if not fluid_solver_settings["compute_reactions"].GetBool(): fluid_solver_settings["compute_reactions"].SetBool(True) warn_msg = '"compute_reactions" is switched off for the fluid-solver, switching it on!' - KM.Logger.PrintWarning("::[NavierStokesAleFluidSolver]::", warn_msg) + KM.Logger.PrintWarning(cls.__name__, warn_msg) else: fluid_solver_settings.AddEmptyValue("compute_reactions").SetBool(True) info_msg = 'Setting "compute_reactions" to true for the fluid-solver' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) @classmethod def _ManipulateMeshMotionSolverSettingsForMeshVelocityComputation(cls, fluid_solver_settings, mesh_motion_solver_settings): @@ -41,7 +39,7 @@ def _ManipulateMeshMotionSolverSettingsForMeshVelocityComputation(cls, fluid_sol if mesh_motion_solver_settings.Has("calculate_mesh_velocity"): if not mesh_motion_solver_settings["calculate_mesh_velocity"].GetBool(): mesh_motion_solver_settings["calculate_mesh_velocity"].SetBool(True) - KM.Logger.PrintWarning("::[NavierStokesAleFluidSolver]::", 'Mesh velocity calculation was deactivated. Switching "calculate_mesh_velocity" on') + KM.Logger.PrintWarning(cls.__name__, 'Mesh velocity calculation was deactivated. Switching "calculate_mesh_velocity" on') else: mesh_motion_solver_settings.AddEmptyValue("calculate_mesh_velocity").SetBool(True) @@ -66,13 +64,13 @@ def _ManipulateMeshMotionSolverSettingsForMeshVelocityComputation(cls, fluid_sol info_msg = '"time_scheme" of the fluid (' + time_scheme_fluid info_msg += ') is different from the\n"time_scheme" used for the ' info_msg += 'calculation of the mesh-velocity (' + time_scheme_mesh_vel + ')' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) else: mesh_vel_calc_settings.AddValue("time_scheme", fluid_solver_settings["time_scheme"]) info_msg = 'setting "time_scheme" of the mesh-solver for the\ncalculation of the ' info_msg += 'mesh-velocity to "' + time_scheme_fluid + '" to be consistent with the\n' info_msg += '"time_scheme" of the fluid' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) if not mesh_vel_calc_settings["time_scheme"].GetString().startswith("bdf"): if mesh_vel_calc_settings.Has("alpha_m"): alpha_mesh_vel = mesh_vel_calc_settings["alpha_m"].GetDouble() @@ -80,13 +78,13 @@ def _ManipulateMeshMotionSolverSettingsForMeshVelocityComputation(cls, fluid_sol info_msg = '"alpha" of the fluid (' + str(alpha_fluid) info_msg += ') is different from the\n"alpha_m" used for the ' info_msg += 'calculation of the mesh-velocity (' + str(alpha_mesh_vel) + ')' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) else: mesh_vel_calc_settings.AddValue("alpha_m", fluid_solver_settings["alpha"]) info_msg = 'setting "alpha_m" of the mesh-solver for the\ncalculation of the ' info_msg += 'mesh-velocity to "' + str(alpha_fluid) + '" to be consistent\nwith ' info_msg += '"alpha" of the fluid' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) elif fluid_solver_type == "fractional_step" or fluid_solver_type == "FractionalStep": # currently fractional step always uses BDF2 @@ -96,17 +94,17 @@ def _ManipulateMeshMotionSolverSettingsForMeshVelocityComputation(cls, fluid_sol info_msg = '"time_scheme" of the fluid (bdf2) ' info_msg += 'is different from the\n"time_scheme" used for the ' info_msg += 'calculation of the mesh-velocity (' + time_scheme_mesh_vel + ')' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) else: mesh_vel_calc_settings.AddEmptyValue("time_scheme").SetString("bdf2") info_msg = 'setting "time_scheme" of the mesh-solver for the\ncalculation of the ' info_msg += 'mesh-velocity to "bdf2" to be consistent with the\n' info_msg += '"time_scheme" of the fluid' - KM.Logger.PrintInfo("::[NavierStokesAleFluidSolver]::", info_msg) + KM.Logger.PrintInfo(cls.__name__, info_msg) if not mesh_vel_calc_settings.Has("time_scheme"): mesh_vel_calc_settings.AddEmptyValue("time_scheme").SetString("UNSPECIFIED") warn_msg = 'unknown "solver_type" of the fluid-solver, therefore ' warn_msg += 'no automatic selection\nof "time_scheme" for the calculation ' warn_msg += 'of the mesh-velocity performed' - KM.Logger.PrintWarning("::[NavierStokesAleFluidSolver]::", warn_msg) + KM.Logger.PrintWarning(cls.__name__, warn_msg) diff --git a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_compressible_solver.py b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_compressible_solver.py index e7aeb9bb79bf..a5f4a65787bd 100644 --- a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_compressible_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_compressible_solver.py @@ -1,4 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.FluidDynamicsApplication as KratosFluid @@ -6,7 +5,6 @@ ## Import base class file from KratosMultiphysics.FluidDynamicsApplication.fluid_solver import FluidSolver -from KratosMultiphysics import python_linear_solver_factory as linear_solver_factory from KratosMultiphysics.FluidDynamicsApplication import check_and_prepare_model_process_fluid def CreateSolver(model, custom_settings): @@ -30,6 +28,7 @@ def GetDefaultSettings(cls): "maximum_iterations": 10, "echo_level": 1, "time_order": 2, + "time_scheme": "bdf2", "compute_reactions": false, "reform_dofs_at_each_step" : true, "relative_tolerance" : 1e-3, @@ -64,20 +63,18 @@ def GetDefaultSettings(cls): def __init__(self, model, custom_settings): self._validate_settings_in_baseclass=True # To be removed eventually + # TODO: DO SOMETHING IN HERE TO REMOVE THE "time_order" FROM THE DEFAULT SETTINGS BUT KEEPING THE BACKWARDS COMPATIBILITY super(NavierStokesCompressibleSolver,self).__init__(model,custom_settings) + self.min_buffer_size = 3 self.element_name = "CompressibleNavierStokes" self.condition_name = "Condition" - self.min_buffer_size = 3 - - ## Construct the linear solver - self.linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) + self.element_integrates_in_time = True ## Set the element replace settings #self._SetCompressibleElementReplaceSettings() - print("Construction of NavierStokesCompressibleSolver finished.") - + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesCompressibleSolver finished.") def AddVariables(self): @@ -108,7 +105,7 @@ def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.PRESSURE) self.main_model_part.AddNodalSolutionStepVariable(KratosFluid.MACH) #for momentum - print("Monolithic compressible fluid solver variables added correctly") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Monolithic compressible fluid solver variables added correctly") def AddDofs(self): KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MOMENTUM_X, KratosMultiphysics.REACTION_X, self.main_model_part) @@ -118,68 +115,21 @@ def AddDofs(self): KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.TOTAL_ENERGY, KratosFluid.REACTION_ENERGY, self.main_model_part) def Initialize(self): - self.computing_model_part = self.GetComputingModelPart() - - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - print("ERROR: _GetAutomaticTimeSteppingUtility out of date") - #self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - # Set the time discretization utility to compute the BDF coefficients - time_order = self.settings["time_order"].GetInt() - if time_order == 2: - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) - else: - raise Exception("Only \"time_order\" equal to 2 is supported. Provided \"time_order\": " + str(time_order)) - - # Creating the solution strategy - self.conv_criteria = KratosMultiphysics.ResidualCriteria(self.settings["relative_tolerance"].GetDouble(), - self.settings["absolute_tolerance"].GetDouble()) - - - #(self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt() - (self.conv_criteria).SetEchoLevel(3) - - domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] - rotation_utility = KratosFluid.CompressibleElementRotationUtility(domain_size,KratosMultiphysics.SLIP) - time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip(rotation_utility) - #time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticScheme() # DOFs (4,5) - - - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(self.computing_model_part, - time_scheme, - self.linear_solver, - self.conv_criteria, - builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - #(self.solver).SetEchoLevel(1) - - - (self.solver).Initialize() - - - # self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) # REMEMBER TO CHECK MY STAB CONSTANTS - - print ("Monolithic compressible solver initialization finished.") + # Construct and set the solution strategy + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) + solution_strategy.Initialize() + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def InitializeSolutionStep(self): (self.time_discretization).ComputeAndSaveBDFCoefficients(self.GetComputingModelPart().ProcessInfo) - (self.solver).InitializeSolutionStep() + self._GetSolutionStrategy().InitializeSolutionStep() def Solve(self): (self.time_discretization).ComputeAndSaveBDFCoefficients(self.GetComputingModelPart().ProcessInfo) - (self.solver).Solve() + self._GetSolutionStrategy().Solve() def PrepareModelPart(self): super(NavierStokesCompressibleSolver,self).PrepareModelPart() @@ -197,24 +147,35 @@ def _ExecuteAfterReading(self): check_and_prepare_model_process_fluid.CheckAndPrepareModelProcess(self.main_model_part, prepare_model_part_settings).Execute() - - #def _SetCompressibleElementReplaceSettings(self): - #domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] - #self.settings.AddEmptyValue("element_replace_settings") - - #if(domain_size == 3): - #self.settings["element_replace_settings"] = KratosMultiphysics.Parameters(""" - #{ - #"element_name":"CompressibleNavierStokes3D4N", - #"condition_name": "SurfaceCondition3D3N" - #} - #""") - #elif(domain_size == 2): - #self.settings["element_replace_settings"] = KratosMultiphysics.Parameters(""" - #{ - #"element_name":"CompressibleNavierStokes2D3N", - #"condition_name": "LineCondition2D2N" - #} - #""") - #else: - #raise Exception("Domain size is not 2 or 3!!") + def _CreateScheme(self): + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + # Cases in which the element manages the time integration + if self.element_integrates_in_time: + # Rotation utility for compressible Navier-Stokes formulations + # A custom rotation util is required as the nodal DOFs differs from the standard incompressible case + rotation_utility = KratosFluid.CompressibleElementRotationUtility( + domain_size, + KratosMultiphysics.SLIP) + # "Fake" scheme for those cases in where the element manages the time integration + # It is required to perform the nodal update once the current time step is solved + scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip(rotation_utility) + # In case the BDF2 scheme is used inside the element, the BDF time discretization utility is required to update the BDF coefficients + if (self.settings["time_scheme"].GetString() == "bdf2"): + time_order = 2 + self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) + else: + err_msg = "Requested elemental time scheme \"" + self.settings["time_scheme"].GetString()+ "\" is not available.\n" + err_msg += "Available options are: \"bdf2\"" + raise Exception(err_msg) + # Cases in which a time scheme manages the time integration + else: + err_msg = "Custom scheme creation is not allowed. Compressible Navier-Stokes elements manage the time integration internally." + raise Exception(err_msg) + return scheme + + def _CreateConvergenceCriterion(self): + convergence_criterion = KratosMultiphysics.ResidualCriteria( + self.settings["relative_tolerance"].GetDouble(), + self.settings["absolute_tolerance"].GetDouble()) + convergence_criterion.SetEchoLevel(self.settings["echo_level"].GetInt()) + return convergence_criterion \ No newline at end of file diff --git a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_embedded_solver.py b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_embedded_solver.py index 4ce6ee0d8419..5042fc6db2a3 100755 --- a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_embedded_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_embedded_solver.py @@ -1,8 +1,5 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics -import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory # Import applications import KratosMultiphysics.FluidDynamicsApplication as KratosCFD @@ -53,6 +50,7 @@ def _SetUpClassicEmbeddedNavierStokes(self, formulation_settings): self.element_name = "EmbeddedNavierStokes" self.condition_name = "NavierStokesWallCondition" self.level_set_type = formulation_settings["level_set_type"].GetString() + self.element_integrates_in_time = True self.element_has_nodal_properties = True self.process_info_data[KratosMultiphysics.DYNAMIC_TAU] = formulation_settings["dynamic_tau"].GetDouble() @@ -74,6 +72,7 @@ def _SetUpEmbeddedSymbolicNavierStokes(self, formulation_settings): self.element_name = "EmbeddedSymbolicNavierStokes" self.condition_name = "NavierStokesWallCondition" self.level_set_type = formulation_settings["level_set_type"].GetString() + self.element_integrates_in_time = True self.element_has_nodal_properties = False self.process_info_data[KratosMultiphysics.DYNAMIC_TAU] = formulation_settings["dynamic_tau"].GetDouble() @@ -94,6 +93,7 @@ def _SetUpClassicEmbeddedAusasNavierStokes(self, formulation_settings): self.element_name = "EmbeddedAusasNavierStokes" self.condition_name = "EmbeddedAusasNavierStokesWallCondition" self.level_set_type = formulation_settings["level_set_type"].GetString() + self.element_integrates_in_time = True self.element_has_nodal_properties = True self.process_info_data[KratosMultiphysics.DYNAMIC_TAU] = formulation_settings["dynamic_tau"].GetDouble() @@ -113,6 +113,7 @@ def _SetUpEmbeddedSymbolicNavierStokesDiscontinuous(self, formulation_settings): self.element_name = "EmbeddedSymbolicNavierStokesDiscontinuous" self.condition_name = "NavierStokesWallCondition" self.level_set_type = formulation_settings["level_set_type"].GetString() + self.element_integrates_in_time = True self.element_has_nodal_properties = False self.process_info_data[KratosMultiphysics.DYNAMIC_TAU] = formulation_settings["dynamic_tau"].GetDouble() @@ -231,8 +232,10 @@ def GetDefaultSettings(cls): "maximum_iterations": 7, "echo_level": 0, "time_order": 2, + "time_scheme": "bdf2", "compute_reactions": false, "reform_dofs_at_each_step": false, + "consider_periodic_conditions": false, "relative_velocity_tolerance": 1e-3, "absolute_velocity_tolerance": 1e-5, "relative_pressure_tolerance": 1e-3, @@ -242,8 +245,8 @@ def GetDefaultSettings(cls): }, "volume_model_part_name" : "volume_model_part", "skin_parts": [""], - "assign_neighbour_elements_to_conditions": false, "no_skin_parts":[""], + "assign_neighbour_elements_to_conditions": false, "time_stepping" : { "automatic_time_step" : true, "CFL_number" : 1, @@ -277,34 +280,28 @@ def ValidateSettings(self): def __init__(self, model, custom_settings): self._validate_settings_in_baseclass=True # To be removed eventually + # TODO: DO SOMETHING IN HERE TO REMOVE THE "time_order" FROM THE DEFAULT SETTINGS BUT KEEPING THE BACKWARDS COMPATIBILITY super(NavierStokesEmbeddedMonolithicSolver,self).__init__(model,custom_settings) self.min_buffer_size = 3 self.embedded_formulation = EmbeddedFormulation(self.settings["formulation"]) self.element_name = self.embedded_formulation.element_name self.condition_name = self.embedded_formulation.condition_name - - ## Set the formulation level set type self.level_set_type = self.embedded_formulation.level_set_type - - ## Set the nodal properties flag + self.element_integrates_in_time = self.embedded_formulation.element_integrates_in_time self.element_has_nodal_properties = self.embedded_formulation.element_has_nodal_properties - ## Construct the linear solver - self.linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - ## Set the distance reading filename # TODO: remove the manual "distance_file_name" set as soon as the problem type one has been tested. if (self.settings["distance_reading_settings"]["import_mode"].GetString() == "from_GiD_file"): self.settings["distance_reading_settings"]["distance_file_name"].SetString(self.settings["model_import_settings"]["input_filename"].GetString()+".post.res") - # If the FM-ALE is required, do a first call to the _get_fm_ale_virtual_model_part + # If the FM-ALE is required, do a first call to the __GetFmAleVirtualModelPart # Note that this will create the virtual model part in the model - self.__fm_ale_is_active = self.settings["fm_ale_settings"]["fm_ale_step_frequency"].GetInt() > 0 - if self.__fm_ale_is_active: - self._get_fm_ale_virtual_model_part() + if self._FmAleIsActive(): + self.__GetFmAleVirtualModelPart() - KratosMultiphysics.Logger.PrintInfo("NavierStokesEmbeddedMonolithicSolver", "Construction of NavierStokesEmbeddedMonolithicSolver finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesEmbeddedMonolithicSolver finished.") def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DENSITY) # TODO: Remove this once the "old" embedded elements get the density from the properties (or once we delete them) @@ -324,18 +321,19 @@ def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosCFD.EMBEDDED_WET_PRESSURE) # Post-process variable (stores the fluid nodes pressure and is set to 0 in the structure ones) self.main_model_part.AddNodalSolutionStepVariable(KratosCFD.EMBEDDED_WET_VELOCITY) # Post-process variable (stores the fluid nodes velocity and is set to 0 in the structure ones) - if self.__fm_ale_is_active: + if self._FmAleIsActive(): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.MESH_DISPLACEMENT) self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.MESH_REACTION) - KratosMultiphysics.Logger.PrintInfo("NavierStokesEmbeddedMonolithicSolver", "Fluid solver variables added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid solver variables added correctly.") def AddDofs(self): super(NavierStokesEmbeddedMonolithicSolver, self).AddDofs() - if self.__fm_ale_is_active: + if self._FmAleIsActive(): KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MESH_DISPLACEMENT_X, KratosMultiphysics.MESH_REACTION_X, self.main_model_part) KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MESH_DISPLACEMENT_Y, KratosMultiphysics.MESH_REACTION_Y, self.main_model_part) KratosMultiphysics.VariableUtils().AddDof(KratosMultiphysics.MESH_DISPLACEMENT_Z, KratosMultiphysics.MESH_REACTION_Z, self.main_model_part) + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "FM-ALE DOFs added correctly.") def PrepareModelPart(self): # Call the base solver PrepareModelPart() @@ -344,77 +342,48 @@ def PrepareModelPart(self): # Set the extra requirements of the embedded formulation if not self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED]: ## Sets the embedded formulation configuration - self._set_embedded_formulation() + self.__SetEmbeddedFormulation() ## Setting the nodal distance - self._set_distance_function() + self.__SetDistanceFunction() def Initialize(self): - computing_model_part = self.GetComputingModelPart() - - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - # Set the time discretization utility to compute the BDF coefficients - time_order = self.settings["time_order"].GetInt() - if time_order == 2: - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) - else: - raise Exception("Only \"time_order\" equal to 2 is supported. Provided \"time_order\": " + str(time_order)) + # If the solver requires an instance of the stabilized embedded_formulation class, set the process info variables + if hasattr(self, 'embedded_formulation'): + self.embedded_formulation.SetProcessInfo(self.GetComputingModelPart()) - # Creating the solution strategy - self.conv_criteria = KratosCFD.VelPrCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) - - time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], # Domain size (2,3) - self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) # DOFs (3,4) - - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(computing_model_part, - time_scheme, - self.linear_solver, - self.conv_criteria, - builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - (self.solver).Initialize() # Initialize the solver. Otherwise the constitutive law is not initializated. + # Construct and initialize the solution strategy + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) + solution_strategy.Initialize() # Set the distance modification process - self._GetDistanceModificationProcess().ExecuteInitialize() + self.GetDistanceModificationProcess().ExecuteInitialize() # For the primitive Ausas formulation, set the find nodal neighbours process # Recall that the Ausas condition requires the nodal neighbours. if (self.settings["formulation"]["element_type"].GetString() == "embedded_ausas_navier_stokes"): - number_of_avg_elems = 10 - number_of_avg_nodes = 10 - self.find_nodal_neighbours_process = KratosMultiphysics.FindNodalNeighboursProcess(self.GetComputingModelPart()) + computing_model_part = self.GetComputingModelPart() + data_communicator = computing_model_part.GetCommunicator().GetDataCommunicator() + self.find_nodal_neighbours_process = KratosMultiphysics.FindGlobalNodalElementalNeighboursProcess( + data_communicator, + computing_model_part) # If required, intialize the FM-ALE utility - if self.__fm_ale_is_active: + if self._FmAleIsActive(): self.fm_ale_step = 1 # Fill the virtual model part geometry. Note that the mesh moving util is created in this first call - self._get_mesh_moving_util().Initialize(self.main_model_part) + self.__GetFmAleUtility().Initialize(self.main_model_part) - KratosMultiphysics.Logger.PrintInfo("NavierStokesEmbeddedMonolithicSolver", "Solver initialization finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def AdvanceInTime(self, current_time): # Call base solver AdvanceInTime to clone the time step and get the new time new_time = super(NavierStokesEmbeddedMonolithicSolver, self).AdvanceInTime(current_time) # Save the current step and time in the virtual model part process info - if self.__fm_ale_is_active: - self._get_fm_ale_virtual_model_part().ProcessInfo[KratosMultiphysics.STEP] += 1 - self._get_fm_ale_virtual_model_part().ProcessInfo[KratosMultiphysics.TIME] = new_time + if self._FmAleIsActive(): + self.__GetFmAleVirtualModelPart().ProcessInfo[KratosMultiphysics.STEP] += 1 + self.__GetFmAleVirtualModelPart().ProcessInfo[KratosMultiphysics.TIME] = new_time return new_time @@ -428,7 +397,7 @@ def InitializeSolutionStep(self): (self.find_nodal_neighbours_process).Execute() # Set the virtual mesh values from the background mesh - self._set_virtual_mesh_values() + self.__SetVirtualMeshValues() # Call the base solver InitializeSolutionStep() super(NavierStokesEmbeddedMonolithicSolver, self).InitializeSolutionStep() @@ -438,11 +407,11 @@ def SolveSolutionStep(self): # Correct the distance field # Note that this is intentionally placed in here (and not in the InitializeSolutionStep() of the solver # It has to be done before each call to the Solve() in case an outer non-linear iteration is performed (FSI) - self._GetDistanceModificationProcess().ExecuteInitializeSolutionStep() + self.GetDistanceModificationProcess().ExecuteInitializeSolutionStep() # Perform the FM-ALE operations # Note that this also sets the EMBEDDED_VELOCITY from the MESH_VELOCITY - self._do_fm_ale_operations() + self.__DoFmAleOperations() # Call the base SolveSolutionStep to solve the embedded CFD problem is_converged = super(NavierStokesEmbeddedMonolithicSolver,self).SolveSolutionStep() @@ -453,7 +422,7 @@ def SolveSolutionStep(self): # Restore the fluid node fixity to its original status # Note that this is intentionally placed in here (and not in the FinalizeSolutionStep() of the solver # It has to be done after each call to the Solve() and the FM-ALE in case an outer non-linear iteration is performed (FSI) - self._GetDistanceModificationProcess().ExecuteFinalizeSolutionStep() + self.GetDistanceModificationProcess().ExecuteFinalizeSolutionStep() return is_converged else: @@ -506,14 +475,7 @@ def _SetNodalProperties(self): KratosMultiphysics.VariableUtils().SetVariable(KratosMultiphysics.DENSITY, rho, self.main_model_part.Nodes) KratosMultiphysics.VariableUtils().SetVariable(KratosMultiphysics.DYNAMIC_VISCOSITY, dyn_viscosity, self.main_model_part.Nodes) - def _set_constitutive_law(self): - ## Construct the constitutive law needed for the embedded element - if(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3): - self.main_model_part.Properties[1][KratosMultiphysics.CONSTITUTIVE_LAW] = KratosCFD.Newtonian3DLaw() - elif(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2): - self.main_model_part.Properties[1][KratosMultiphysics.CONSTITUTIVE_LAW] = KratosCFD.Newtonian2DLaw() - - def _set_embedded_formulation(self): + def __SetEmbeddedFormulation(self): # Set the SLIP elemental flag if (self.settings["formulation"]["is_slip"].GetBool()): KratosMultiphysics.VariableUtils().SetFlag(KratosMultiphysics.SLIP, True, self.GetComputingModelPart().Elements) @@ -524,7 +486,7 @@ def _set_embedded_formulation(self): # Save the formulation settings in the ProcessInfo self.embedded_formulation.SetProcessInfo(self.main_model_part) - def _set_distance_function(self): + def __SetDistanceFunction(self): ## Set the nodal distance function if (self.settings["distance_reading_settings"]["import_mode"].GetString() == "from_GiD_file"): DistanceUtility = read_distance_from_file.DistanceImportUtility(self.main_model_part, self.settings["distance_reading_settings"]) @@ -536,7 +498,7 @@ def _set_distance_function(self): distance_value = node.GetSolutionStepValue(KratosMultiphysics.DISTANCE) node.SetSolutionStepValue(KratosMultiphysics.DISTANCE, -distance_value) - def _GetDistanceModificationProcess(self): + def GetDistanceModificationProcess(self): if not hasattr(self, '_distance_modification_process'): self._distance_modification_process = self.__CreateDistanceModificationProcess() return self._distance_modification_process @@ -550,29 +512,32 @@ def __CreateDistanceModificationProcess(self): distance_modification_settings["model_part_name"].SetString(aux_full_volume_part_name) return KratosCFD.DistanceModificationProcess(self.model, distance_modification_settings) - def _get_fm_ale_structure_model_part(self): + def _FmAleIsActive(self): + return self.settings["fm_ale_settings"]["fm_ale_step_frequency"].GetInt() > 0 + + def __GetFmAleStructureModelPart(self): structure_model_part_name = self.settings["fm_ale_settings"]["fm_ale_solver_settings"]["structure_model_part_name"].GetString() if self.model.HasModelPart(structure_model_part_name): return self.model.GetModelPart(structure_model_part_name) else: raise Exception("Structure model part {0} not found in model.\n It is expected to be added in your custom analysis stage file.".format(structure_model_part_name)) - def _get_fm_ale_virtual_model_part(self): + def __GetFmAleVirtualModelPart(self): if not hasattr(self, '_virtual_model_part'): - self._virtual_model_part = self._create_fm_ale_virtual_model_part() + self._virtual_model_part = self.__CreateFmAleVirtualModelPart() return self._virtual_model_part - def _create_fm_ale_virtual_model_part(self): + def __CreateFmAleVirtualModelPart(self): virtual_model_part_name = self.settings["fm_ale_settings"]["fm_ale_solver_settings"]["virtual_model_part_name"].GetString() virtual_model_part = self.model.CreateModelPart(virtual_model_part_name) return virtual_model_part - def _get_mesh_moving_util(self): + def __GetFmAleUtility(self): if not hasattr (self, '_mesh_moving_util'): - self._mesh_moving_util = self._create_mesh_moving_util() + self._mesh_moving_util = self.__CreateFmAleUtility() return self._mesh_moving_util - def _create_mesh_moving_util(self): + def __CreateFmAleUtility(self): if have_mesh_moving: mesh_movement = self.settings["fm_ale_settings"]["mesh_movement"].GetString() if (mesh_movement == "implicit"): @@ -590,8 +555,8 @@ def _create_mesh_moving_util(self): else: raise Exception("MeshMovingApplication is required to construct the FM-ALE utility (ExplicitFixedMeshALEUtilities)") - def _is_fm_ale_step(self): - if self.__fm_ale_is_active: + def __IsFmAleStep(self): + if self._FmAleIsActive(): if (self.fm_ale_step == self.settings["fm_ale_settings"]["fm_ale_step_frequency"].GetInt()): return True else: @@ -600,33 +565,33 @@ def _is_fm_ale_step(self): return False def __UpdateFMALEStepCounter(self): - if self.__fm_ale_is_active: - if (self._is_fm_ale_step()): + if self._FmAleIsActive(): + if (self.__IsFmAleStep()): # Reset the FM-ALE steps counter self.fm_ale_step = 1 else: # Update the FM-ALE steps counter self.fm_ale_step += 1 - def _set_virtual_mesh_values(self): - if self._is_fm_ale_step(): + def __SetVirtualMeshValues(self): + if self.__IsFmAleStep(): # Fill the virtual model part variable values: VELOCITY (n,nn), PRESSURE (n,nn) - self._get_mesh_moving_util().SetVirtualMeshValuesFromOriginMesh() + self.__GetFmAleUtility().SetVirtualMeshValuesFromOriginMesh() - def _do_fm_ale_operations(self): - if self._is_fm_ale_step(): + def __DoFmAleOperations(self): + if self.__IsFmAleStep(): # Solve the mesh problem delta_time = self.main_model_part.ProcessInfo[KratosMultiphysics.DELTA_TIME] - self._get_mesh_moving_util().ComputeMeshMovement(delta_time) + self.__GetFmAleUtility().ComputeMeshMovement(delta_time) # Project the obtained MESH_VELOCITY and historical VELOCITY and PRESSURE values to the origin mesh buffer_size = self.main_model_part.GetBufferSize() domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] if (domain_size == 2): - self._get_mesh_moving_util().ProjectVirtualValues2D(self.main_model_part, buffer_size) + self.__GetFmAleUtility().ProjectVirtualValues2D(self.main_model_part, buffer_size) else: - self._get_mesh_moving_util().ProjectVirtualValues3D(self.main_model_part, buffer_size) + self.__GetFmAleUtility().ProjectVirtualValues3D(self.main_model_part, buffer_size) # If FM-ALE is performed, use the MESH_VELOCITY as EMBEDDED_VELOCITY KratosMultiphysics.VariableUtils().CopyModelPartNodalVarToNonHistoricalVar( @@ -637,6 +602,6 @@ def _do_fm_ale_operations(self): 0) def __UndoFMALEOperations(self): - if self._is_fm_ale_step(): + if self.__IsFmAleStep(): # Undo the FM-ALE virtual mesh movement - self._get_mesh_moving_util().UndoMeshMovement() + self.__GetFmAleUtility().UndoMeshMovement() diff --git a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_fractionalstep.py b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_fractionalstep.py index 5f7db2b98141..d2160987c342 100755 --- a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_fractionalstep.py +++ b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_fractionalstep.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory @@ -99,17 +97,16 @@ def __init__(self, model, custom_settings): self.element_name = custom_settings["formulation"]["element_type"].GetString() self.condition_name = custom_settings["formulation"]["condition_type"].GetString() - self.min_buffer_size = 3 self.element_has_nodal_properties = True - ## Construct the linear solvers - self.pressure_linear_solver = linear_solver_factory.ConstructSolver(self.settings["pressure_linear_solver_settings"]) - self.velocity_linear_solver = linear_solver_factory.ConstructSolver(self.settings["velocity_linear_solver_settings"]) + self.min_buffer_size = 3 self.compute_reactions = self.settings["compute_reactions"].GetBool() - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStep", "Construction of NavierStokesSolverFractionalStep solver finished.") + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesSolverFractionalStep solver finished.") def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DENSITY) @@ -137,69 +134,98 @@ def AddVariables(self): if self.settings["consider_periodic_conditions"].GetBool() == True: self.main_model_part.AddNodalSolutionStepVariable(KratosCFD.PATCH_INDEX) - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStep", "Fluid solver variables added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid solver variables added correctly.") def Initialize(self): - self.computing_model_part = self.GetComputingModelPart() - - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - #TODO: next part would be much cleaner if we passed directly the parameters to the c++ - if self.settings["consider_periodic_conditions"] == True: - self.solver_settings = KratosCFD.FractionalStepSettingsPeriodic(self.computing_model_part, - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.settings.GetInt(), - self.settings["use_slip_conditions"].GetBool(), - self.settings["move_mesh_flag"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - KratosCFD.PATCH_INDEX) - - else: - self.solver_settings = KratosCFD.FractionalStepSettings(self.computing_model_part, - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.settings["time_order"].GetInt(), - self.settings["use_slip_conditions"].GetBool(), - self.settings["move_mesh_flag"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool()) - - self.solver_settings.SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.solver_settings.SetStrategy(KratosCFD.StrategyLabel.Velocity, - self.velocity_linear_solver, - self.settings["velocity_tolerance"].GetDouble(), - self.settings["maximum_velocity_iterations"].GetInt()) + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) + solution_strategy.Initialize() - self.solver_settings.SetStrategy(KratosCFD.StrategyLabel.Pressure, - self.pressure_linear_solver, - self.settings["pressure_tolerance"].GetDouble(), - self.settings["maximum_pressure_iterations"].GetInt()) - - - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.solver = KratosCFD.FSStrategy(self.computing_model_part, - self.solver_settings, - self.settings["predictor_corrector"].GetBool(), - KratosCFD.PATCH_INDEX) - else: - self.solver = KratosCFD.FSStrategy(self.computing_model_part, - self.solver_settings, - self.settings["predictor_corrector"].GetBool()) - - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) - - (self.solver).Initialize() - - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverFractionalStep", "Solver initialization finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def SolveSolutionStep(self): if self._TimeBufferIsInitialized(): is_converged = super(NavierStokesSolverFractionalStep,self).SolveSolutionStep() if self.compute_reactions: - self.solver.CalculateReactions() + self._GetSolutionStrategy().CalculateReactions() return is_converged else: return True + + def _CreateScheme(self): + pass + + def _CreateLinearSolver(self): + # Create the pressure linear solver + pressure_linear_solver_configuration = self.settings["pressure_linear_solver_settings"] + pressure_linear_solver = linear_solver_factory.ConstructSolver(pressure_linear_solver_configuration) + # Create the velocity linear solver + velocity_linear_solver_configuration = self.settings["velocity_linear_solver_settings"] + velocity_linear_solver = linear_solver_factory.ConstructSolver(velocity_linear_solver_configuration) + # Return a tuple containing both linear solvers + return (pressure_linear_solver, velocity_linear_solver) + + def _CreateConvergenceCriterion(self): + pass + + def _CreateBuilderAndSolver(self): + pass + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + domain_size = computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + + # Create the pressure and velocity linear solvers + # Note that linear_solvers is a tuple. The first item is the pressure + # linear solver. The second item is the velocity linear solver. + linear_solvers = self._GetLinearSolver() + + # Create the fractional step settings instance + # TODO: next part would be much cleaner if we passed directly the parameters to the c++ + if self.settings["consider_periodic_conditions"].GetBool(): + fractional_step_settings = KratosCFD.FractionalStepSettingsPeriodic( + computing_model_part, + domain_size, + self.settings["time_order"].GetInt(), + self.settings["use_slip_conditions"].GetBool(), + self.settings["move_mesh_flag"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + KratosCFD.PATCH_INDEX) + else: + fractional_step_settings = KratosCFD.FractionalStepSettings( + computing_model_part, + domain_size, + self.settings["time_order"].GetInt(), + self.settings["use_slip_conditions"].GetBool(), + self.settings["move_mesh_flag"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool()) + + # Set the strategy echo level + fractional_step_settings.SetEchoLevel(self.settings["echo_level"].GetInt()) + + # Set the velocity and pressure fractional step strategy settings + fractional_step_settings.SetStrategy(KratosCFD.StrategyLabel.Pressure, + linear_solvers[0], + self.settings["pressure_tolerance"].GetDouble(), + self.settings["maximum_pressure_iterations"].GetInt()) + + fractional_step_settings.SetStrategy(KratosCFD.StrategyLabel.Velocity, + linear_solvers[1], + self.settings["velocity_tolerance"].GetDouble(), + self.settings["maximum_velocity_iterations"].GetInt()) + + # Create the fractional step strategy + if self.settings["consider_periodic_conditions"].GetBool() == True: + solution_strategy = KratosCFD.FSStrategy( + computing_model_part, + fractional_step_settings, + self.settings["predictor_corrector"].GetBool(), + KratosCFD.PATCH_INDEX) + else: + solution_strategy = KratosCFD.FSStrategy( + computing_model_part, + fractional_step_settings, + self.settings["predictor_corrector"].GetBool()) + + return solution_strategy diff --git a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_vmsmonolithic.py b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_vmsmonolithic.py index a27849d3dfa0..dfcfa2292d4e 100755 --- a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_vmsmonolithic.py +++ b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_solver_vmsmonolithic.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics @@ -9,8 +7,7 @@ # Import base class file from KratosMultiphysics.FluidDynamicsApplication.fluid_solver import FluidSolver -import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory - +# Import turbulence model solver from KratosMultiphysics.FluidDynamicsApplication.turbulence_model_solver import CreateTurbulenceModel class StabilizedFormulation(object): @@ -275,17 +272,13 @@ def __init__(self, model, custom_settings): msg += "Accepted values are \"bossak\", \"bdf2\" or \"steady\".\n" raise Exception(msg) - ## Construct the linear solver - self.linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - ## Construct the turbulence model solver if not self.settings["turbulence_model_solver_settings"].IsEquivalentTo(KratosMultiphysics.Parameters("{}")): self._turbulence_model_solver = CreateTurbulenceModel(model, self.settings["turbulence_model_solver_settings"]) self.condition_name = self._turbulence_model_solver.GetFluidVelocityPressureConditionName() - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithic", "Using " + self.condition_name + " as wall condition") - - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithic", "Construction of NavierStokesSolverMonolithic finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Using " + self.condition_name + " as wall condition") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesSolverMonolithic finished.") def AddVariables(self): ## Add base class variables @@ -314,10 +307,11 @@ def AddVariables(self): self._turbulence_model_solver.fluid_model_part = self.main_model_part self._turbulence_model_solver.AddVariables() + # Adding variables required for the periodic conditions if self.settings["consider_periodic_conditions"].GetBool() == True: self.main_model_part.AddNodalSolutionStepVariable(KratosCFD.PATCH_INDEX) - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithic", "Fluid solver variables added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid solver variables added correctly.") def AddDofs(self): super(NavierStokesSolverMonolithic, self).AddDofs() @@ -334,111 +328,20 @@ def PrepareModelPart(self): self._turbulence_model_solver.PrepareModelPart() def Initialize(self): + # If the solver requires an instance of the stabilized formulation class, set the process info variables + if hasattr(self, 'formulation'): + self.formulation.SetProcessInfo(self.GetComputingModelPart()) - self.computing_model_part = self.GetComputingModelPart() - - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - # Creating the solution strategy - - if self.settings["time_scheme"].GetString() == "steady": - self.conv_criteria = KratosMultiphysics.ResidualCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble()) - else: - self.conv_criteria = KratosCFD.VelPrCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) - - # Creating the time integration scheme - if (self.element_integrates_in_time): - # "Fake" scheme for those cases in where the element manages the time integration - # It is required to perform the nodal update once the current time step is solved - self.time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip( - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) - # In case the BDF2 scheme is used inside the element, the BDF time discretization utility is required to update the BDF coefficients - if (self.settings["time_scheme"].GetString() == "bdf2"): - time_order = 2 - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) - else: - err_msg = "Requested elemental time scheme \"" + self.settings["time_scheme"].GetString()+ "\" is not available.\n" - err_msg += "Available options are: \"bdf2\"" - raise Exception(err_msg) - else: - if not hasattr(self, "_turbulence_model_solver"): - # Bossak time integration scheme - if self.settings["time_scheme"].GetString() == "bossak": - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.time_scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( - self.settings["alpha"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - KratosCFD.PATCH_INDEX) - else: - self.time_scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( - self.settings["alpha"].GetDouble(), - self.settings["move_mesh_strategy"].GetInt(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - # BDF2 time integration scheme - elif self.settings["time_scheme"].GetString() == "bdf2": - self.time_scheme = KratosCFD.GearScheme() - # Time scheme for steady state fluid solver - elif self.settings["time_scheme"].GetString() == "steady": - self.time_scheme = KratosCFD.ResidualBasedSimpleSteadyScheme( - self.settings["velocity_relaxation"].GetDouble(), - self.settings["pressure_relaxation"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - else: - err_msg = "Requested time scheme " + self.settings["time_scheme"].GetString() + " is not available.\n" - err_msg += "Available options are: \"bossak\", \"bdf2\" and \"steady\"" - raise Exception(err_msg) - else: - self._turbulence_model_solver.Initialize() - if self.settings["time_scheme"].GetString() == "bossak": - self.time_scheme = KratosCFD.ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent( - self.settings["alpha"].GetDouble(), - self.settings["move_mesh_strategy"].GetInt(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.settings["turbulence_model_solver_settings"]["velocity_pressure_relaxation_factor"].GetDouble(), - self._turbulence_model_solver.GetTurbulenceSolvingProcess()) - # Time scheme for steady state fluid solver - elif self.settings["time_scheme"].GetString() == "steady": - self.time_scheme = KratosCFD.ResidualBasedSimpleSteadyScheme( - self.settings["velocity_relaxation"].GetDouble(), - self.settings["pressure_relaxation"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self._turbulence_model_solver.GetTurbulenceSolvingProcess()) - if self.settings["consider_periodic_conditions"].GetBool(): - builder_and_solver = KratosCFD.ResidualBasedBlockBuilderAndSolverPeriodic(self.linear_solver, - KratosCFD.PATCH_INDEX) - else: - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(self.computing_model_part, - self.time_scheme, - self.linear_solver, - self.conv_criteria, - builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.formulation.SetProcessInfo(self.computing_model_part) - - (self.solver).Initialize() + # Construct and initialize the solution strategy + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) + solution_strategy.Initialize() + # If there is turbulence modelling, set the new solution strategy as parent strategy if hasattr(self, "_turbulence_model_solver"): - self._turbulence_model_solver.SetParentSolvingStrategy(self.solver) + self._turbulence_model_solver.SetParentSolvingStrategy(solution_strategy) - KratosMultiphysics.Logger.PrintInfo("NavierStokesSolverMonolithic", "Solver initialization finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def InitializeSolutionStep(self): if self._TimeBufferIsInitialized(): @@ -446,7 +349,7 @@ def InitializeSolutionStep(self): if hasattr(self, 'time_discretization'): (self.time_discretization).ComputeAndSaveBDFCoefficients(self.GetComputingModelPart().ProcessInfo) # Perform the solver InitializeSolutionStep - (self.solver).InitializeSolutionStep() + self._GetSolutionStrategy().InitializeSolutionStep() # Perform the turbulence modelling InitializeSolutionStep if hasattr(self, "_turbulence_model_solver"): self._turbulence_model_solver.InitializeSolutionStep() diff --git a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_two_fluids_solver.py b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_two_fluids_solver.py index 7945e8596ac3..aeff79f90c3f 100644 --- a/applications/FluidDynamicsApplication/python_scripts/navier_stokes_two_fluids_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/navier_stokes_two_fluids_solver.py @@ -1,5 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.kratos_utilities as KratosUtilities @@ -14,8 +12,6 @@ from KratosMultiphysics.FluidDynamicsApplication.fluid_solver import FluidSolver from KratosMultiphysics.FluidDynamicsApplication.read_distance_from_file import DistanceImportUtility -import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory - def CreateSolver(model, custom_settings): return NavierStokesTwoFluidsSolver(model, custom_settings) @@ -44,8 +40,10 @@ def GetDefaultSettings(cls): "maximum_iterations": 7, "echo_level": 0, "time_order": 2, + "time_scheme": "bdf2", "compute_reactions": false, "reform_dofs_at_each_step": false, + "consider_periodic_conditions": false, "relative_velocity_tolerance": 1e-3, "absolute_velocity_tolerance": 1e-5, "relative_pressure_tolerance": 1e-3, @@ -55,7 +53,7 @@ def GetDefaultSettings(cls): }, "volume_model_part_name" : "volume_model_part", "skin_parts": [""], - "assign_neighbour_elements_to_conditions": false, + "assign_neighbour_elements_to_conditions": true, "no_skin_parts":[""], "time_stepping" : { "automatic_time_step" : true, @@ -78,22 +76,29 @@ def GetDefaultSettings(cls): def __init__(self, model, custom_settings): self._validate_settings_in_baseclass=True # To be removed eventually + + # TODO: DO SOMETHING IN HERE TO REMOVE THE "time_order" FROM THE DEFAULT SETTINGS BUT KEEPING THE BACKWARDS COMPATIBILITY + super(NavierStokesTwoFluidsSolver,self).__init__(model,custom_settings) self.element_name = "TwoFluidNavierStokes" self.condition_name = "NavierStokesWallCondition" + self.element_integrates_in_time = True self.element_has_nodal_properties = True self.min_buffer_size = 3 self._bfecc_convection = self.settings["bfecc_convection"].GetBool() + dynamic_tau = self.settings["formulation"]["dynamic_tau"].GetDouble() + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, dynamic_tau) + ## Set the distance reading filename # TODO: remove the manual "distance_file_name" set as soon as the problem type one has been tested. if (self.settings["distance_reading_settings"]["import_mode"].GetString() == "from_GiD_file"): self.settings["distance_reading_settings"]["distance_file_name"].SetString(self.settings["model_import_settings"]["input_filename"].GetString()+".post.res") - KratosMultiphysics.Logger.PrintInfo("NavierStokesTwoFluidsSolver", "Construction of NavierStokesTwoFluidsSolver finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of NavierStokesTwoFluidsSolver finished.") def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DENSITY) @@ -114,76 +119,38 @@ def AddVariables(self): self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DISTANCE) # Distance function nodal values self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DISTANCE_GRADIENT) # Distance gradient nodal values - KratosMultiphysics.Logger.PrintInfo("NavierStokesTwoFluidsSolver", "Fluid solver variables added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Fluid solver variables added correctly.") def PrepareModelPart(self): # Initialize the level-set function if not self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED]: ## Setting the nodal distance - self._set_distance_function() + self.__SetDistanceFunction() # Call the base solver PrepareModelPart() super(NavierStokesTwoFluidsSolver, self).PrepareModelPart() def Initialize(self): - self.computing_model_part = self.GetComputingModelPart() - - ## Construct the linear solver - self.linear_solver = linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - - KratosMultiphysics.NormalCalculationUtils().CalculateOnSimplex(self.computing_model_part, self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - - self.neighbour_search = KratosMultiphysics.FindNodalNeighboursProcess(self.computing_model_part) - (self.neighbour_search).Execute() - - self.accelerationLimitationUtility = KratosCFD.AccelerationLimitationUtilities( self.computing_model_part, 5.0 ) - - # If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - # Set the time discretization utility to compute the BDF coefficients - time_order = self.settings["time_order"].GetInt() - if time_order == 2: - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) - else: - raise Exception("Only \"time_order\" equal to 2 is supported. Provided \"time_order\": " + str(time_order)) - - # Creating the solution strategy - self.conv_criteria = KratosCFD.VelPrCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) + computing_model_part = self.GetComputingModelPart() - self.level_set_convection_process = self._set_level_set_convection_process() + # Calculate boundary normals + KratosMultiphysics.NormalCalculationUtils().CalculateOnSimplex( + computing_model_part, + computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - self.variational_distance_process = self._set_variational_distance_process() + # Elemental neighbours calculation + data_communicator = computing_model_part.GetCommunicator().GetDataCommunicator() + neighbour_search = KratosMultiphysics.FindGlobalNodalElementalNeighboursProcess( + data_communicator, + computing_model_part) + neighbour_search.Execute() - time_scheme = KratosMultiphysics.ResidualBasedIncrementalUpdateStaticSchemeSlip(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], # Domain size (2,3) - self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) # DOFs (3,4) + # Set and initialize the solution strategy + solution_strategy = self._GetSolutionStrategy() + solution_strategy.SetEchoLevel(self.settings["echo_level"].GetInt()) + solution_strategy.Initialize() - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(self.computing_model_part, - time_scheme, - self.linear_solver, - self.conv_criteria, - builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - (self.solver).Initialize() # Initialize the solver. Otherwise the constitutive law is not initializated. - (self.solver).Check() - - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["formulation"]["dynamic_tau"].GetDouble()) - - KratosMultiphysics.Logger.PrintInfo("NavierStokesTwoFluidsSolver", "Solver initialization finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def InitializeSolutionStep(self): if self._TimeBufferIsInitialized(): @@ -192,31 +159,30 @@ def InitializeSolutionStep(self): # Perform the level-set convection according to the previous step velocity if self._bfecc_convection: - (self.level_set_convection_process).BFECCconvect( + self._GetLevelSetConvectionProcess().BFECCconvect( self.main_model_part, KratosMultiphysics.DISTANCE, KratosMultiphysics.VELOCITY, self.settings["bfecc_number_substeps"].GetInt()) else: - (self.level_set_convection_process).Execute() + self._GetLevelSetConvectionProcess().Execute() # Recompute the distance field according to the new level-set position - (self.variational_distance_process).Execute() + self._GetVariationalDistanceProcess().Execute() # Update the DENSITY and DYNAMIC_VISCOSITY values according to the new level-set self._SetNodalProperties() # Initialize the solver current step - (self.solver).InitializeSolutionStep() + self._GetSolutionStrategy().InitializeSolutionStep() def FinalizeSolutionStep(self): if self._TimeBufferIsInitialized(): - (self.solver).FinalizeSolutionStep() - (self.accelerationLimitationUtility).Execute() + self._GetSolutionStrategy().FinalizeSolutionStep() + self._GetAccelerationLimitationUtility().Execute() # TODO: Remove this method as soon as the subproperties are available def _SetPhysicalProperties(self): - import os warn_msg = '\nThe materials import mechanism used in the two fluids solver is DEPRECATED!\n' warn_msg += 'It will be removed to use the base fluid_solver.py one as soon as the subproperties are available.\n' KratosMultiphysics.Logger.PrintWarning('\n\x1b[1;31mDEPRECATION-WARNING\x1b[0m', warn_msg) @@ -242,8 +208,7 @@ def _SetPhysicalProperties(self): aux_material_settings = KratosMultiphysics.Parameters("""{"Parameters": {"materials_filename": ""}} """) aux_material_settings["Parameters"]["materials_filename"].SetString(aux_materials_filename) KratosMultiphysics.ReadMaterialsUtility(aux_material_settings, self.model) - - os.remove(aux_materials_filename) + KratosUtilities.DeleteFileIfExisting(aux_materials_filename) materials_imported = True else: @@ -284,7 +249,7 @@ def _SetNodalProperties(self): node.SetSolutionStepValue(KratosMultiphysics.DENSITY, rho_2) node.SetSolutionStepValue(KratosMultiphysics.DYNAMIC_VISCOSITY, mu_2) - def _set_distance_function(self): + def __SetDistanceFunction(self): ## Set the nodal distance function if (self.settings["distance_reading_settings"]["import_mode"].GetString() == "from_GiD_file"): DistanceUtility = DistanceImportUtility(self.main_model_part, self.settings["distance_reading_settings"]) @@ -292,45 +257,73 @@ def _set_distance_function(self): elif (self.settings["distance_reading_settings"]["import_mode"].GetString() == "from_mdpa"): KratosMultiphysics.Logger.PrintInfo("Navier Stokes Embedded Solver","Distance function taken from the .mdpa input file.") - def _set_level_set_convection_process(self): + def _GetAccelerationLimitationUtility(self): + if not hasattr(self, '_acceleration_limitation_utility'): + self._acceleration_limitation_utility = self.__CreateAccelerationLimitationUtility() + return self._acceleration_limitation_utility + + def _GetLevelSetConvectionProcess(self): + if not hasattr(self, '_level_set_convection_process'): + self._level_set_convection_process = self._CreateLevelSetConvectionProcess() + return self._level_set_convection_process + + def _GetVariationalDistanceProcess(self): + if not hasattr(self, '_variational_distance_process'): + self._variational_distance_process = self._CreateVariationalDistanceProcess() + return self._variational_distance_process + + def __CreateAccelerationLimitationUtility(self): + maximum_multiple_of_g_acceleration_allowed = 5.0 + acceleration_limitation_utility = KratosCFD.AccelerationLimitationUtilities( + self.GetComputingModelPart(), + maximum_multiple_of_g_acceleration_allowed) + + return acceleration_limitation_utility + + def _CreateLevelSetConvectionProcess(self): # Construct the level set convection process + domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + computing_model_part = self.GetComputingModelPart() if self._bfecc_convection: if have_conv_diff: - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: - locator = KratosMultiphysics.BinBasedFastPointLocator2D(self.main_model_part).UpdateSearchDatabase() + if domain_size == 2: + locator = KratosMultiphysics.BinBasedFastPointLocator2D(computing_model_part).UpdateSearchDatabase() level_set_convection_process = KratosConvDiff.BFECCConvection2D(locator) else: - locator = KratosMultiphysics.BinBasedFastPointLocator3D(self.main_model_part).UpdateSearchDatabase() + locator = KratosMultiphysics.BinBasedFastPointLocator3D(computing_model_part).UpdateSearchDatabase() level_set_convection_process = KratosConvDiff.BFECCConvection3D(locator) else: raise Exception("The BFECC level set convection requires the Kratos ConvectionDiffusionApplication compilation.") else: - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: + linear_solver = self._GetLinearSolver() + if domain_size == 2: level_set_convection_process = KratosMultiphysics.LevelSetConvectionProcess2D( KratosMultiphysics.DISTANCE, - self.main_model_part, - self.linear_solver) + computing_model_part, + linear_solver) else: level_set_convection_process = KratosMultiphysics.LevelSetConvectionProcess3D( KratosMultiphysics.DISTANCE, - self.main_model_part, - self.linear_solver) + computing_model_part, + linear_solver) return level_set_convection_process - def _set_variational_distance_process(self): + def _CreateVariationalDistanceProcess(self): # Construct the variational distance calculation process maximum_iterations = 2 #TODO: Make this user-definable + linear_solver = self._GetLinearSolver() + computing_model_part = self.GetComputingModelPart() if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: variational_distance_process = KratosMultiphysics.VariationalDistanceCalculationProcess2D( - self.main_model_part, - self.linear_solver, + computing_model_part, + linear_solver, maximum_iterations, KratosMultiphysics.VariationalDistanceCalculationProcess2D.CALCULATE_EXACT_DISTANCES_TO_PLANE) else: variational_distance_process = KratosMultiphysics.VariationalDistanceCalculationProcess3D( - self.main_model_part, - self.linear_solver, + computing_model_part, + linear_solver, maximum_iterations, KratosMultiphysics.VariationalDistanceCalculationProcess3D.CALCULATE_EXACT_DISTANCES_TO_PLANE) diff --git a/applications/FluidDynamicsApplication/python_scripts/steady_navier_stokes_solver_vmsmonolithic.py b/applications/FluidDynamicsApplication/python_scripts/steady_navier_stokes_solver_vmsmonolithic.py deleted file mode 100644 index 198066f2f8bf..000000000000 --- a/applications/FluidDynamicsApplication/python_scripts/steady_navier_stokes_solver_vmsmonolithic.py +++ /dev/null @@ -1,87 +0,0 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 -# importing the Kratos Library -import KratosMultiphysics -import KratosMultiphysics.FluidDynamicsApplication as KratosCFD - -## Import base class file -from KratosMultiphysics.FluidDynamicsApplication.navier_stokes_solver_vmsmonolithic import NavierStokesSolverMonolithic - -def CreateSolver(model, custom_settings): - return SteadyNavierStokesSolver_VMSMonolithic(model, custom_settings) - -class SteadyNavierStokesSolver_VMSMonolithic(NavierStokesSolverMonolithic): - - def __init__(self, model, custom_settings): - - # parse and strip parameters that do not exist in base class. we need to remove - # extra parameters so base class doesn't throw an error. alternatively a single solver script - # could be used and the scheme type could be passed in json parameters. - self.velocity_relaxation_factor = custom_settings["velocity_relaxation_factor"].GetDouble() - self.pressure_relaxation_factor = custom_settings["pressure_relaxation_factor"].GetDouble() - base_settings = custom_settings - base_settings.RemoveValue("velocity_relaxation_factor") - base_settings.RemoveValue("pressure_relaxation_factor") - - # call base class constructor with remaining parameters - super().__init__(model, base_settings) - - if self.settings["consider_periodic_conditions"].GetBool() == True: - raise ValueError("consider_periodic_conditions not supported yet.") - - if self.settings["move_mesh_strategy"].GetInt() != 0: - raise ValueError("move_mesh_strategy not supported yet.") - - def Initialize(self): - - self.computing_model_part = self.GetComputingModelPart() - - MoveMeshFlag = False - - # TODO: TURBULENCE MODELS ARE NOT ADDED YET - #~ # Turbulence model - #~ if self.use_spalart_allmaras: - #~ self.activate_spalart_allmaras() - - # creating the solution strategy - self.conv_criteria = KratosCFD.VelPrCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.time_scheme = KratosCFD.ResidualBasedSimpleSteadyScheme(self.velocity_relaxation_factor, self.pressure_relaxation_factor, self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - - # TODO: TURBULENCE MODELS ARE NOT ADDED YET - #~ if self.turbulence_model is None: - #~ if self.periodic == True: - #~ self.time_scheme = ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent\ - #~ (self.alpha, self.move_mesh_strategy, self.domain_size, PATCH_INDEX) - #~ else: - #~ self.time_scheme = ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent\ - #~ (self.alpha, self.move_mesh_strategy, self.domain_size) - #~ else: - #~ self.time_scheme = ResidualBasedPredictorCorrectorVelocityBossakSchemeTurbulent\ - #~ (self.alpha, self.move_mesh_strategy, self.domain_size, self.turbulence_model) - - builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(self.linear_solver) - - - self.solver = KratosMultiphysics.ResidualBasedNewtonRaphsonStrategy(self.main_model_part, - self.time_scheme, - self.linear_solver, - self.conv_criteria, - builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - if self.settings["stabilization"].Has("dynamic_tau"): - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["stabilization"]["dynamic_tau"].GetDouble()) - if self.settings["stabilization"].Has("oss_switch"): - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["stabilization"]["oss_switch"].GetInt()) - - print ("Monolithic solver initialization finished.") diff --git a/applications/FluidDynamicsApplication/python_scripts/stokes_solver.py b/applications/FluidDynamicsApplication/python_scripts/stokes_solver.py index 8c5d418dc02d..d4a2347f5ec2 100644 --- a/applications/FluidDynamicsApplication/python_scripts/stokes_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/stokes_solver.py @@ -1,4 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 # importing the Kratos Library import KratosMultiphysics as kratoscore import KratosMultiphysics.FluidDynamicsApplication as cfd diff --git a/applications/FluidDynamicsApplication/python_scripts/stokes_solver_monolithic.py b/applications/FluidDynamicsApplication/python_scripts/stokes_solver_monolithic.py index 0b5688d59ff6..9ab0d3b747aa 100644 --- a/applications/FluidDynamicsApplication/python_scripts/stokes_solver_monolithic.py +++ b/applications/FluidDynamicsApplication/python_scripts/stokes_solver_monolithic.py @@ -1,5 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # importing the Kratos Library import KratosMultiphysics as kratoscore import KratosMultiphysics.python_linear_solver_factory as linear_solver_factory diff --git a/applications/FluidDynamicsApplication/python_scripts/trilinos_adjoint_vmsmonolithic_solver.py b/applications/FluidDynamicsApplication/python_scripts/trilinos_adjoint_vmsmonolithic_solver.py index fbbc17feb564..a7416a25dd23 100644 --- a/applications/FluidDynamicsApplication/python_scripts/trilinos_adjoint_vmsmonolithic_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/trilinos_adjoint_vmsmonolithic_solver.py @@ -1,4 +1,3 @@ -from __future__ import print_function, absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 ## Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.mpi as KratosMPI @@ -6,7 +5,6 @@ # Import applications import KratosMultiphysics.TrilinosApplication as TrilinosApplication from KratosMultiphysics.TrilinosApplication import trilinos_linear_solver_factory -import KratosMultiphysics.FluidDynamicsApplication as FluidDynamicsApplication from KratosMultiphysics.FluidDynamicsApplication.adjoint_vmsmonolithic_solver import AdjointVMSMonolithicSolver from KratosMultiphysics.mpi.distributed_import_model_part_utility import DistributedImportModelPartUtility @@ -37,7 +35,7 @@ def GetDefaultSettings(cls): "materials_filename": "" }, "linear_solver_settings" : { - "solver_type" : "multi_level" + "solver_type" : "amgcl" }, "volume_model_part_name" : "volume_model_part", "skin_parts": [""], @@ -53,7 +51,9 @@ def GetDefaultSettings(cls): "time_stepping": { "automatic_time_step" : false, "time_step" : -0.1 - } + }, + "consider_periodic_conditions": false, + "assign_neighbour_elements_to_conditions": false }""") default_settings.AddMissingParameters(super(AdjointVMSMonolithicMPISolver, cls).GetDefaultSettings()) @@ -71,8 +71,8 @@ def __init__(self, model, custom_settings): self.min_buffer_size = 2 self.element_has_nodal_properties = True - # construct the linear solver - self.trilinos_linear_solver = trilinos_linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of AdjointVMSMonolithicMPISolver finished.") @@ -91,71 +91,73 @@ def ImportModelPart(self): ## Execute the Metis partitioning and reading self.distributed_model_part_importer.ExecutePartitioningAndReading() - KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "TrilinosNavierStokesSolverMonolithic","MPI model reading finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "MPI model reading finished.") def PrepareModelPart(self): super(self.__class__,self).PrepareModelPart() ## Construct the MPI communicators self.distributed_model_part_importer.CreateCommunicators() - def AddDofs(self): - super(self.__class__, self).AddDofs() - - KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "DOFs for the AdjointVMSMonolithicMPISolver added correctly in all processors.") - - def Initialize(self): - - ## Construct the communicator - self.EpetraCommunicator = TrilinosApplication.CreateCommunicator() - - ## Get the computing model part - self.computing_model_part = self.GetComputingModelPart() - - domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] - if self.settings["response_function_settings"]["response_type"].GetString() == "drag": - if (domain_size == 2): - self.response_function = FluidDynamicsApplication.DragResponseFunction2D(self.settings["response_function_settings"]["custom_settings"], self.main_model_part) - elif (domain_size == 3): - self.response_function = FluidDynamicsApplication.DragResponseFunction3D(self.settings["response_function_settings"]["custom_settings"], self.main_model_part) - else: - raise Exception("Invalid DOMAIN_SIZE: " + str(domain_size)) + def _GetEpetraCommunicator(self): + if not hasattr(self, '_epetra_communicator'): + self._epetra_communicator = TrilinosApplication.CreateCommunicator() + return self._epetra_communicator + + def _CreateScheme(self): + response_function = self.GetResponseFunction() + scheme_type = self.settings["scheme_settings"]["scheme_type"].GetString() + if scheme_type == "bossak": + scheme = TrilinosApplication.TrilinosResidualBasedAdjointBossakScheme( + self.settings["scheme_settings"], + response_function) + elif scheme_type == "steady": + scheme = TrilinosApplication.TrilinosResidualBasedAdjointSteadyScheme(response_function) else: - raise Exception("invalid response_type: " + self.settings["response_function_settings"]["response_type"].GetString()) + raise Exception("Invalid scheme_type: " + scheme_type) + return scheme - self.sensitivity_builder = KratosMultiphysics.SensitivityBuilder(self.settings["sensitivity_settings"], self.main_model_part, self.response_function) - - if self.settings["scheme_settings"]["scheme_type"].GetString() == "bossak": - self.time_scheme = TrilinosApplication.TrilinosResidualBasedAdjointBossakScheme(self.settings["scheme_settings"], self.response_function) - elif self.settings["scheme_settings"]["scheme_type"].GetString() == "steady": - self.time_scheme = TrilinosApplication.TrilinosResidualBasedAdjointSteadyScheme(self.response_function) - else: - raise Exception("invalid scheme_type: " + self.settings["scheme_settings"]["scheme_type"].GetString()) + def _CreateLinearSolver(self): + linear_solver_configuration = self.settings["linear_solver_settings"] + return trilinos_linear_solver_factory.ConstructSolver(linear_solver_configuration) - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3: + def _CreateBuilderAndSolver(self): + # Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + if domain_size == 3: guess_row_size = 20*4 - elif self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: + else: guess_row_size = 10*3 - - self.builder_and_solver = TrilinosApplication.TrilinosBlockBuilderAndSolver(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver) - - self.solver = TrilinosApplication.TrilinosLinearStrategy(self.main_model_part, - self.time_scheme, - self.trilinos_linear_solver, - self.builder_and_solver, - False, - False, - False, - False) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - (self.solver).Initialize() - (self.response_function).Initialize() - (self.sensitivity_builder).Initialize() - - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) - - KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"Monolithic MPI solver initialization finished.") + # Construct the Trilinos builder and solver + trilinos_linear_solver = self._GetLinearSolver() + epetra_communicator = self._GetEpetraCommunicator() + if self.settings["consider_periodic_conditions"].GetBool(): + builder_and_solver = TrilinosApplication.TrilinosBlockBuilderAndSolverPeriodic( + epetra_communicator, + guess_row_size, + trilinos_linear_solver, + KratosFluid.PATCH_INDEX) + else: + builder_and_solver = TrilinosApplication.TrilinosBlockBuilderAndSolver( + epetra_communicator, + guess_row_size, + trilinos_linear_solver) + return builder_and_solver + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + builder_and_solver = self._GetBuilderAndSolver() + calculate_reaction_flag = False + reform_dof_set_at_each_step = False + calculate_norm_dx_flag = False + move_mesh_flag = False + return TrilinosApplication.TrilinosLinearStrategy( + computing_model_part, + time_scheme, + linear_solver, + builder_and_solver, + calculate_reaction_flag, + reform_dof_set_at_each_step, + calculate_norm_dx_flag, + move_mesh_flag) diff --git a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_embedded_solver.py b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_embedded_solver.py index 78477e84509e..537e82276566 100644 --- a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_embedded_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_embedded_solver.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.mpi as KratosMPI # MPI-python interface @@ -19,7 +17,6 @@ def CreateSolver(model, custom_settings): class NavierStokesMPIEmbeddedMonolithicSolver(navier_stokes_embedded_solver.NavierStokesEmbeddedMonolithicSolver): - @classmethod def GetDefaultSettings(cls): @@ -44,20 +41,15 @@ def GetDefaultSettings(cls): "echo_level": 0, "consider_periodic_conditions": false, "time_order": 2, + "time_scheme": "bdf2", "compute_reactions": false, "reform_dofs_at_each_step": false, "relative_velocity_tolerance": 1e-3, "absolute_velocity_tolerance": 1e-5, "relative_pressure_tolerance": 1e-3, "absolute_pressure_tolerance": 1e-5, - "linear_solver_settings" : { - "solver_type" : "multi_level", - "max_iteration" : 200, - "tolerance" : 1e-6, - "max_levels" : 3, - "symmetric" : false, - "reform_preconditioner_at_each_step" : true, - "scaling" : true + "linear_solver_settings": { + "solver_type": "amgcl" }, "volume_model_part_name" : "volume_model_part", "skin_parts": [""], @@ -68,7 +60,6 @@ def GetDefaultSettings(cls): "minimum_delta_time" : 1e-4, "maximum_delta_time" : 0.01 }, - "periodic": "periodic", "move_mesh_flag": false, "formulation": { "element_type": "embedded_element_from_defaults", @@ -87,29 +78,24 @@ def GetDefaultSettings(cls): def __init__(self, model, custom_settings): self._validate_settings_in_baseclass=True # To be removed eventually # Note: deliberately calling the constructor of the base python solver (the parent of my parent) + # TODO: ONCE THE FM-ALE WORKS IN MPI, IT WOULD BE POSSIBLE TO ONLY CALL THE BASE CLASS CONSTRUCTOR super(navier_stokes_embedded_solver.NavierStokesEmbeddedMonolithicSolver, self).__init__(model,custom_settings) self.min_buffer_size = 3 self.embedded_formulation = navier_stokes_embedded_solver.EmbeddedFormulation(self.settings["formulation"]) self.element_name = self.embedded_formulation.element_name self.condition_name = self.embedded_formulation.condition_name - - ## Set the formulation level set type self.level_set_type = self.embedded_formulation.level_set_type - - ## Set the nodal properties flag + self.element_integrates_in_time = self.embedded_formulation.element_integrates_in_time self.element_has_nodal_properties = self.embedded_formulation.element_has_nodal_properties - ## Construct the linear solver - self.trilinos_linear_solver = trilinos_linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPIEmbeddedMonolithicSolver","Construction of NavierStokesMPIEmbeddedMonolithicSolver finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"Construction of NavierStokesMPIEmbeddedMonolithicSolver finished.") def AddVariables(self): super(NavierStokesMPIEmbeddedMonolithicSolver, self).AddVariables() self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.PARTITION_INDEX) - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPIEmbeddedMonolithicSolver","Variables for the Trilinos monolithic embedded fluid solver added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"Variables for the Trilinos monolithic embedded fluid solver added correctly.") def ImportModelPart(self): ## Construct the Distributed import model part utility @@ -118,89 +104,94 @@ def ImportModelPart(self): self.distributed_model_part_importer.ImportModelPart() ## Sets DENSITY, VISCOSITY and SOUND_VELOCITY - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPIEmbeddedMonolithicSolver","MPI model reading finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"MPI model reading finished.") def PrepareModelPart(self): super(NavierStokesMPIEmbeddedMonolithicSolver,self).PrepareModelPart() ## Construct MPI the communicators self.distributed_model_part_importer.CreateCommunicators() - def AddDofs(self): - super(NavierStokesMPIEmbeddedMonolithicSolver, self).AddDofs() - - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPIEmbeddedMonolithicSolver","DOFs for the VMS Trilinos fluid solver added correctly.") - - def Initialize(self): - ## Construct the communicator - self.EpetraCommunicator = KratosTrilinos.CreateCommunicator() - - ## Get the computing model part - self.computing_model_part = self.GetComputingModelPart() - - ## If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - # Set the time discretization utility to compute the BDF coefficients - time_order = self.settings["time_order"].GetInt() - if time_order == 2: - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) + def _GetEpetraCommunicator(self): + if not hasattr(self, '_epetra_communicator'): + self._epetra_communicator = KratosTrilinos.CreateCommunicator() + return self._epetra_communicator + + def _CreateScheme(self): + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + # Cases in which the element manages the time integration + if self.element_integrates_in_time: + # "Fake" scheme for those cases in where the element manages the time integration + # It is required to perform the nodal update once the current time step is solved + scheme = KratosTrilinos.TrilinosResidualBasedIncrementalUpdateStaticSchemeSlip( + domain_size, + domain_size + 1) + # In case the BDF2 scheme is used inside the element, the BDF time discretization utility is required to update the BDF coefficients + if (self.settings["time_scheme"].GetString() == "bdf2"): + time_order = 2 + self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) + else: + err_msg = "Requested elemental time scheme \"" + self.settings["time_scheme"].GetString()+ "\" is not available.\n" + err_msg += "Available options are: \"bdf2\"" + raise Exception(err_msg) + # Cases in which a time scheme manages the time integration else: - raise Exception("Only \"time_order\" equal to 2 is supported. Provided \"time_order\": " + str(time_order)) - - ## Creating the Trilinos convergence criteria - self.conv_criteria = KratosTrilinos.TrilinosUPCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) - - ## Creating the Trilinos incremental update time scheme (the time integration is defined within the embedded element) - self.time_scheme = KratosTrilinos.TrilinosResidualBasedIncrementalUpdateStaticSchemeSlip(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], # Domain size (2,3) - self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) # DOFs (3,4) - - ## Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3: + err_msg = "Custom scheme creation is not allowed. Embedded Navier-Stokes elements manage the time integration internally." + raise Exception(err_msg) + return scheme + + def _CreateLinearSolver(self): + linear_solver_configuration = self.settings["linear_solver_settings"] + return trilinos_linear_solver_factory.ConstructSolver(linear_solver_configuration) + + def _CreateConvergenceCriterion(self): + convergence_criterion = KratosTrilinos.TrilinosUPCriteria( + self.settings["relative_velocity_tolerance"].GetDouble(), + self.settings["absolute_velocity_tolerance"].GetDouble(), + self.settings["relative_pressure_tolerance"].GetDouble(), + self.settings["absolute_pressure_tolerance"].GetDouble()) + convergence_criterion.SetEchoLevel(self.settings["echo_level"].GetInt()) + return convergence_criterion + + def _CreateBuilderAndSolver(self): + # Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + if domain_size == 3: guess_row_size = 20*4 - elif self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: + else: guess_row_size = 10*3 - - ## Construct the Trilinos builder and solver - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolverPeriodic(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver, - KratosFluid.PATCH_INDEX) + # Construct the Trilinos builder and solver + trilinos_linear_solver = self._GetLinearSolver() + epetra_communicator = self._GetEpetraCommunicator() + if self.settings["consider_periodic_conditions"].GetBool(): + builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolverPeriodic( + epetra_communicator, + guess_row_size, + trilinos_linear_solver, + KratosFluid.PATCH_INDEX) else: - self.builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolver(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver) - - ## Construct the Trilinos Newton-Raphson strategy - self.solver = KratosTrilinos.TrilinosNewtonRaphsonStrategy(self.main_model_part, - self.time_scheme, - self.trilinos_linear_solver, - self.conv_criteria, - self.builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - (self.solver).Initialize() - - # Set the distance modification process - self.__GetDistanceModificationProcess().ExecuteInitialize() - - # For the primitive Ausas formulation, set the find nodal neighbours process - # Recall that the Ausas condition requires the nodal neighbouts. - if (self.settings["formulation"]["element_type"].GetString() == "embedded_ausas_navier_stokes"): - number_of_avg_elems = 10 - number_of_avg_nodes = 10 - self.find_nodal_neighbours_process = KratosMultiphysics.FindNodalNeighboursProcess(self.GetComputingModelPart(), - number_of_avg_elems, - number_of_avg_nodes) - - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPIEmbeddedMonolithicSolver","Solver initialization finished.") \ No newline at end of file + builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolver( + epetra_communicator, + guess_row_size, + trilinos_linear_solver) + return builder_and_solver + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + convergence_criterion = self._GetConvergenceCriterion() + builder_and_solver = self._GetBuilderAndSolver() + return KratosTrilinos.TrilinosNewtonRaphsonStrategy( + computing_model_part, + time_scheme, + linear_solver, + convergence_criterion, + builder_and_solver, + self.settings["maximum_iterations"].GetInt(), + self.settings["compute_reactions"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + self.settings["move_mesh_flag"].GetBool()) + + @classmethod + def _FmAleIsActive(self): + return False \ No newline at end of file diff --git a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_fractionalstep.py b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_fractionalstep.py index 14804c58d4a4..b1f6d4d4c352 100755 --- a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_fractionalstep.py +++ b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_fractionalstep.py @@ -1,12 +1,9 @@ -from __future__ import absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.mpi as KratosMPI # MPI-python interface # Import applications import KratosMultiphysics.TrilinosApplication as KratosTrilinos # MPI solvers -import KratosMultiphysics.FluidDynamicsApplication as KratosFluid # Fluid dynamics application from KratosMultiphysics.FluidDynamicsApplication import TrilinosExtension as TrilinosFluid from KratosMultiphysics.TrilinosApplication import trilinos_linear_solver_factory @@ -47,22 +44,10 @@ def GetDefaultSettings(cls): "compute_reactions": false, "reform_dofs_at_each_step": false, "pressure_linear_solver_settings": { - "solver_type" : "multi_level", - "max_iteration" : 200, - "tolerance" : 1e-6, - "symmetric" : true, - "scaling" : true, - "reform_preconditioner_at_each_step" : true, - "verbosity" : 0 + "solver_type": "amgcl" }, "velocity_linear_solver_settings": { - "solver_type" : "multi_level", - "max_iteration" : 200, - "tolerance" : 1e-6, - "symmetric" : false, - "scaling" : true, - "reform_preconditioner_at_each_step" : true, - "verbosity" : 0 + "solver_type": "amgcl" }, "volume_model_part_name" : "volume_model_part", "skin_parts":[""], @@ -90,13 +75,12 @@ def __init__(self, model, custom_settings): self.min_buffer_size = 3 self.element_has_nodal_properties = True - ## Construct the linear solvers - self.pressure_linear_solver = trilinos_linear_solver_factory.ConstructSolver(self.settings["pressure_linear_solver_settings"]) - self.velocity_linear_solver = trilinos_linear_solver_factory.ConstructSolver(self.settings["velocity_linear_solver_settings"]) - self.compute_reactions = self.settings["compute_reactions"].GetBool() - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverFractionalStep","Construction of TrilinosNavierStokesSolverFractionalStep solver finished.") + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) + + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"Construction of TrilinosNavierStokesSolverFractionalStep solver finished.") def AddVariables(self): @@ -105,9 +89,8 @@ def AddVariables(self): ## Add specific MPI variables self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.PARTITION_INDEX) - self.main_model_part.AddNodalSolutionStepVariable(KratosFluid.PATCH_INDEX) - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverFractionalStep","variables for the trilinos fractional step solver added correctly") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"variables for the trilinos fractional step solver added correctly") def ImportModelPart(self): @@ -116,7 +99,7 @@ def ImportModelPart(self): ## Execute the Metis partitioning and reading self.distributed_model_part_importer.ImportModelPart() - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverFractionalStep","MPI model reading finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"MPI model reading finished.") def PrepareModelPart(self): super(TrilinosNavierStokesSolverFractionalStep,self).PrepareModelPart() @@ -128,63 +111,99 @@ def AddDofs(self): ## Base class DOFs addition super(TrilinosNavierStokesSolverFractionalStep, self).AddDofs() - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverFractionalStep","DOFs for the VMS Trilinos fluid solver added correctly in all processors.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"DOFs for the VMS Trilinos fluid solver added correctly in all processors.") def Initialize(self): - ## Construct the communicator - self.EpetraComm = KratosTrilinos.CreateCommunicator() - - ## Get the computing model part - self.computing_model_part = self.GetComputingModelPart() - - ## If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - #TODO: next part would be much cleaner if we passed directly the parameters to the c++ - if self.settings["consider_periodic_conditions"] == True: - self.solver_settings = TrilinosFluid.TrilinosFractionalStepSettingsPeriodic(self.EpetraComm, - self.computing_model_part, - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.settings["time_order"].GetInt(), - self.settings["use_slip_conditions"].GetBool(), - self.settings["move_mesh_flag"].GetBool(), - self.settings["reform_dofs_at_each_step]"].GetBool(), - KratosFluid.PATCH_INDEX) - - else: - self.solver_settings = TrilinosFluid.TrilinosFractionalStepSettings(self.EpetraComm, - self.computing_model_part, - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.settings["time_order"].GetInt(), - self.settings["use_slip_conditions"].GetBool(), - self.settings["move_mesh_flag"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool()) - - self.solver_settings.SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.solver_settings.SetStrategy(TrilinosFluid.TrilinosStrategyLabel.Velocity, - self.velocity_linear_solver, - self.settings["velocity_tolerance"].GetDouble(), - self.settings["maximum_velocity_iterations"].GetInt()) - - self.solver_settings.SetStrategy(TrilinosFluid.TrilinosStrategyLabel.Pressure, - self.pressure_linear_solver, - self.settings["pressure_tolerance"].GetDouble(), - self.settings["maximum_pressure_iterations"].GetInt()) - - self.solver = TrilinosFluid.TrilinosFSStrategy(self.computing_model_part, - self.solver_settings, - self.settings["predictor_corrector"].GetBool(), - KratosFluid.PATCH_INDEX) + ## Base class solver intiialization + super(TrilinosNavierStokesSolverFractionalStep, self).Initialize() - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["dynamic_tau"].GetDouble()) - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.OSS_SWITCH, self.settings["oss_switch"].GetInt()) - - (self.solver).Initialize() - - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverFractionalStep","Initialization TrilinosNavierStokesSolverFractionalStep finished") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") def Finalize(self): - self.solver.Clear() + self._GetSolutionStrategy().Clear() + + def _GetEpetraCommunicator(self): + if not hasattr(self, '_epetra_communicator'): + self._epetra_communicator = KratosTrilinos.CreateCommunicator() + return self._epetra_communicator + + def _CreateScheme(self): + pass + + def _CreateLinearSolver(self): + # Create the pressure linear solver + pressure_linear_solver_configuration = self.settings["pressure_linear_solver_settings"] + pressure_linear_solver = trilinos_linear_solver_factory.ConstructSolver(pressure_linear_solver_configuration) + # Create the velocity linear solver + velocity_linear_solver_configuration = self.settings["velocity_linear_solver_settings"] + velocity_linear_solver = trilinos_linear_solver_factory.ConstructSolver(velocity_linear_solver_configuration) + # Return a tuple containing both linear solvers + return (pressure_linear_solver, velocity_linear_solver) + + def _CreateConvergenceCriterion(self): + pass + + def _CreateBuilderAndSolver(self): + pass + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + epetra_communicator = self._GetEpetraCommunicator() + domain_size = computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + + # Create the pressure and velocity linear solvers + # Note that linear_solvers is a tuple. The first item is the pressure + # linear solver. The second item is the velocity linear solver. + linear_solvers = self._GetLinearSolver() + + # Create the fractional step settings instance + # TODO: next part would be much cleaner if we passed directly the parameters to the c++ + if self.settings["consider_periodic_conditions"].GetBool(): + fractional_step_settings = TrilinosFluid.TrilinosFractionalStepSettingsPeriodic( + epetra_communicator, + computing_model_part, + domain_size, + self.settings["time_order"].GetInt(), + self.settings["use_slip_conditions"].GetBool(), + self.settings["move_mesh_flag"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + KratosCFD.PATCH_INDEX) + else: + fractional_step_settings = TrilinosFluid.TrilinosFractionalStepSettings( + epetra_communicator, + computing_model_part, + domain_size, + self.settings["time_order"].GetInt(), + self.settings["use_slip_conditions"].GetBool(), + self.settings["move_mesh_flag"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool()) + + # Set the strategy echo level + fractional_step_settings.SetEchoLevel(self.settings["echo_level"].GetInt()) + + # Set the velocity and pressure fractional step strategy settings + fractional_step_settings.SetStrategy(TrilinosFluid.TrilinosStrategyLabel.Pressure, + linear_solvers[0], + self.settings["pressure_tolerance"].GetDouble(), + self.settings["maximum_pressure_iterations"].GetInt()) + + fractional_step_settings.SetStrategy(TrilinosFluid.TrilinosStrategyLabel.Velocity, + linear_solvers[1], + self.settings["velocity_tolerance"].GetDouble(), + self.settings["maximum_velocity_iterations"].GetInt()) + + # Create the fractional step strategy + if self.settings["consider_periodic_conditions"].GetBool(): + solution_strategy = TrilinosFluid.TrilinosFSStrategy( + computing_model_part, + fractional_step_settings, + self.settings["predictor_corrector"].GetBool(), + KratosCFD.PATCH_INDEX) + else: + solution_strategy = TrilinosFluid.TrilinosFSStrategy( + computing_model_part, + fractional_step_settings, + self.settings["predictor_corrector"].GetBool()) + + return solution_strategy diff --git a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_vmsmonolithic.py b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_vmsmonolithic.py index 21bf2bbbca80..fb7bd4be090e 100755 --- a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_vmsmonolithic.py +++ b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_solver_vmsmonolithic.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.mpi as KratosMPI # MPI-python interface @@ -41,21 +39,14 @@ def GetDefaultSettings(cls): "maximum_iterations": 10, "echo_level": 0, "consider_periodic_conditions": false, - "time_order": 2, "compute_reactions": false, "reform_dofs_at_each_step": false, "relative_velocity_tolerance": 1e-5, "absolute_velocity_tolerance": 1e-7, "relative_pressure_tolerance": 1e-5, "absolute_pressure_tolerance": 1e-7, - "linear_solver_settings" : { - "solver_type" : "multi_level", - "max_iteration" : 200, - "tolerance" : 1e-8, - "max_levels" : 3, - "symmetric" : false, - "reform_preconditioner_at_each_step" : true, - "scaling" : true + "linear_solver_settings": { + "solver_type": "amgcl" }, "volume_model_part_name" : "volume_model_part", "skin_parts": [""], @@ -104,17 +95,13 @@ def __init__(self, model, custom_settings): msg += "Accepted values are \"bossak\", \"bdf2\" or \"steady\".\n" raise Exception(msg) - ## Construct the linear solver - self.trilinos_linear_solver = trilinos_linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) - ## Construct the turbulence model solver if not self.settings["turbulence_model_solver_settings"].IsEquivalentTo(KratosMultiphysics.Parameters("{}")): self._turbulence_model_solver = CreateTurbulenceModel(model, self.settings["turbulence_model_solver_settings"]) self.condition_name = self._turbulence_model_solver.GetFluidVelocityPressureConditionName() - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverMonolithic", "Using " + self.condition_name + " as wall condition") - - KratosMultiphysics.Logger.Print("Construction of TrilinosNavierStokesSolverMonolithic finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Using " + self.condition_name + " as wall condition") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Construction of TrilinosNavierStokesSolverMonolithic finished.") def AddVariables(self): ## Add variables from the base class @@ -131,146 +118,151 @@ def ImportModelPart(self): ## Execute the Metis partitioning and reading self.distributed_model_part_importer.ImportModelPart() - KratosMultiphysics.Logger.PrintInfo("TrilinosNavierStokesSolverMonolithic","MPI model reading finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"MPI model reading finished.") def PrepareModelPart(self): - if not self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED]: - self._SetPhysicalProperties() - - super(navier_stokes_solver_vmsmonolithic.NavierStokesSolverMonolithic, self).PrepareModelPart() + # Call the base solver to do the PrepareModelPart + # Note that his also calls the PrepareModelPart of the turbulence model + super(TrilinosNavierStokesSolverMonolithic, self).PrepareModelPart() + # Create the MPI communicators self.distributed_model_part_importer.CreateCommunicators() + def Initialize(self): + # If there is turbulence modelling, set the Epetra communicator in the turbulence solver if hasattr(self, "_turbulence_model_solver"): - self._turbulence_model_solver.PrepareModelPart() + self._turbulence_model_solver.SetCommunicator(self._GetEpetraCommunicator()) - def AddDofs(self): - ## Base class DOFs addition - super(TrilinosNavierStokesSolverMonolithic, self).AddDofs() + # Call the base Initialize() method to create and initialize the strategy + super(TrilinosNavierStokesSolverMonolithic, self).Initialize() - KratosMultiphysics.Logger.Print("DOFs for the VMS Trilinos fluid solver added correctly in all processors.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, "Solver initialization finished.") + def Finalize(self): + self._GetSolutionStrategy().Clear() - def Initialize(self): - ## Construct the communicator - self.EpetraCommunicator = KratosTrilinos.CreateCommunicator() if hasattr(self, "_turbulence_model_solver"): - self._turbulence_model_solver.SetCommunicator(self.EpetraCommunicator) - - ## Get the computing model part - self.computing_model_part = self.GetComputingModelPart() - - ## If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - ## Creating the Trilinos convergence criteria - if (self.settings["time_scheme"].GetString() == "bossak"): - self.conv_criteria = KratosTrilinos.TrilinosUPCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - elif (self.settings["time_scheme"].GetString() == "steady"): - self.conv_criteria = KratosTrilinos.TrilinosResidualCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble()) + self._turbulence_model_solver.Finalize() - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) + def _GetEpetraCommunicator(self): + if not hasattr(self, '_epetra_communicator'): + self._epetra_communicator = KratosTrilinos.CreateCommunicator() + return self._epetra_communicator - ## Creating the Trilinos time scheme - if (self.element_integrates_in_time): + def _CreateScheme(self): + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + # Cases in which the element manages the time integration + if self.element_integrates_in_time: # "Fake" scheme for those cases in where the element manages the time integration # It is required to perform the nodal update once the current time step is solved - self.time_scheme = KratosTrilinos.TrilinosResidualBasedIncrementalUpdateStaticSchemeSlip( - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) + scheme = KratosTrilinos.TrilinosResidualBasedIncrementalUpdateStaticSchemeSlip( + domain_size, + domain_size + 1) # In case the BDF2 scheme is used inside the element, set the time discretization utility to compute the BDF coefficients if (self.settings["time_scheme"].GetString() == "bdf2"): - time_order = self.settings["time_order"].GetInt() - if time_order == 2: - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) - else: - raise Exception("Only \"time_order\" equal to 2 is supported. Provided \"time_order\": " + str(time_order)) + time_order = 2 + self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) else: err_msg = "Requested elemental time scheme " + self.settings["time_scheme"].GetString() + " is not available.\n" err_msg += "Available options are: \"bdf2\"" raise Exception(err_msg) + # Cases in which a time scheme manages the time integration else: if not hasattr(self, "_turbulence_model_solver"): + # Bossak time integration scheme if self.settings["time_scheme"].GetString() == "bossak": # TODO: Can we remove this periodic check, Is the PATCH_INDEX used in this scheme? if self.settings["consider_periodic_conditions"].GetBool() == True: - self.time_scheme = TrilinosFluid.TrilinosPredictorCorrectorVelocityBossakSchemeTurbulent( + scheme = TrilinosFluid.TrilinosPredictorCorrectorVelocityBossakSchemeTurbulent( self.settings["alpha"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], + domain_size, KratosCFD.PATCH_INDEX) else: - self.time_scheme = TrilinosFluid.TrilinosPredictorCorrectorVelocityBossakSchemeTurbulent( + scheme = TrilinosFluid.TrilinosPredictorCorrectorVelocityBossakSchemeTurbulent( self.settings["alpha"].GetDouble(), self.settings["move_mesh_strategy"].GetInt(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) + domain_size) + # BDF2 time integration scheme + elif self.settings["time_scheme"].GetString() == "bdf2": + scheme = TrilinosFluid.TrilinosGearScheme() + # Time scheme for steady state fluid solver elif self.settings["time_scheme"].GetString() == "steady": - self.time_scheme = TrilinosFluid.TrilinosResidualBasedSimpleSteadyScheme( + scheme = TrilinosFluid.TrilinosResidualBasedSimpleSteadyScheme( self.settings["velocity_relaxation"].GetDouble(), self.settings["pressure_relaxation"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) + domain_size) else: self._turbulence_model_solver.Initialize() if self.settings["time_scheme"].GetString() == "bossak": - self.time_scheme = TrilinosFluid.TrilinosPredictorCorrectorVelocityBossakSchemeTurbulent( + scheme = TrilinosFluid.TrilinosPredictorCorrectorVelocityBossakSchemeTurbulent( self.settings["alpha"].GetDouble(), self.settings["move_mesh_strategy"].GetInt(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], + domain_size, self.settings["turbulence_model_solver_settings"]["velocity_pressure_relaxation_factor"].GetDouble(), self._turbulence_model_solver.GetTurbulenceSolvingProcess()) # Time scheme for steady state fluid solver elif self.settings["time_scheme"].GetString() == "steady": - self.time_scheme = TrilinosFluid.TrilinosResidualBasedSimpleSteadyScheme( + scheme = TrilinosFluid.TrilinosResidualBasedSimpleSteadyScheme( self.settings["velocity_relaxation"].GetDouble(), self.settings["pressure_relaxation"].GetDouble(), - self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], + domain_size, self._turbulence_model_solver.GetTurbulenceSolvingProcess()) + return scheme + + def _CreateLinearSolver(self): + linear_solver_configuration = self.settings["linear_solver_settings"] + return trilinos_linear_solver_factory.ConstructSolver(linear_solver_configuration) - ## Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3: + def _CreateConvergenceCriterion(self): + if self.settings["time_scheme"].GetString() == "steady": + convergence_criterion = KratosTrilinos.TrilinosResidualCriteria( + self.settings["relative_velocity_tolerance"].GetDouble(), + self.settings["absolute_velocity_tolerance"].GetDouble()) + else: + convergence_criterion = KratosTrilinos.TrilinosUPCriteria( + self.settings["relative_velocity_tolerance"].GetDouble(), + self.settings["absolute_velocity_tolerance"].GetDouble(), + self.settings["relative_pressure_tolerance"].GetDouble(), + self.settings["absolute_pressure_tolerance"].GetDouble()) + convergence_criterion.SetEchoLevel(self.settings["echo_level"].GetInt()) + return convergence_criterion + + def _CreateBuilderAndSolver(self): + # Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + if domain_size == 3: guess_row_size = 20*4 - elif self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: + else: guess_row_size = 10*3 - - ## Construct the Trilinos builder and solver - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolverPeriodic(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver, - KratosCFD.PATCH_INDEX) + # Construct the Trilinos builder and solver + trilinos_linear_solver = self._GetLinearSolver() + epetra_communicator = self._GetEpetraCommunicator() + if self.settings["consider_periodic_conditions"].GetBool(): + builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolverPeriodic( + epetra_communicator, + guess_row_size, + trilinos_linear_solver, + KratosFluid.PATCH_INDEX) else: - self.builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolver(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver) - - ## Construct the Trilinos Newton-Raphson strategy - self.solver = KratosTrilinos.TrilinosNewtonRaphsonStrategy(self.main_model_part, - self.time_scheme, - self.trilinos_linear_solver, - self.conv_criteria, - self.builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - - self.formulation.SetProcessInfo(self.computing_model_part) - - (self.solver).Initialize() - - if hasattr(self, "_turbulence_model_solver"): - self._turbulence_model_solver.SetParentSolvingStrategy(self.solver) - - KratosMultiphysics.Logger.Print("Monolithic MPI solver initialization finished.") - - def Finalize(self): - self.solver.Clear() - - if hasattr(self, "_turbulence_model_solver"): - self._turbulence_model_solver.Finalize() + builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolver( + epetra_communicator, + guess_row_size, + trilinos_linear_solver) + return builder_and_solver + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + convergence_criterion = self._GetConvergenceCriterion() + builder_and_solver = self._GetBuilderAndSolver() + return KratosTrilinos.TrilinosNewtonRaphsonStrategy( + computing_model_part, + time_scheme, + linear_solver, + convergence_criterion, + builder_and_solver, + self.settings["maximum_iterations"].GetInt(), + self.settings["compute_reactions"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + self.settings["move_mesh_flag"].GetBool()) \ No newline at end of file diff --git a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_two_fluids_solver.py b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_two_fluids_solver.py index 4157aa1173d3..92c4f47cb4b8 100755 --- a/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_two_fluids_solver.py +++ b/applications/FluidDynamicsApplication/python_scripts/trilinos_navier_stokes_two_fluids_solver.py @@ -1,5 +1,3 @@ -from __future__ import absolute_import, division # makes KratosMultiphysics backward compatible with python 2.6 and 2.7 - # Importing the Kratos Library import KratosMultiphysics import KratosMultiphysics.mpi as KratosMPI # MPI-python interface @@ -39,22 +37,17 @@ def GetDefaultSettings(cls): }, "maximum_iterations": 7, "echo_level": 0, - "consider_periodic_conditions" : false, "time_order": 2, + "time_scheme": "bdf2", "compute_reactions": false, "reform_dofs_at_each_step": false, + "consider_periodic_conditions": false, "relative_velocity_tolerance": 1e-3, "absolute_velocity_tolerance": 1e-5, "relative_pressure_tolerance": 1e-3, "absolute_pressure_tolerance": 1e-5, "linear_solver_settings": { - "solver_type": "AmgclMPISolver", - "tolerance": 1.0e-6, - "max_iteration": 10, - "scaling": false, - "verbosity": 0, - "preconditioner_type": "None", - "krylov_type": "cg" + "solver_type": "amgcl" }, "volume_model_part_name" : "volume_model_part", "skin_parts": [""], @@ -68,7 +61,9 @@ def GetDefaultSettings(cls): "move_mesh_flag": false, "formulation": { "dynamic_tau": 1.0 - } + }, + "bfecc_convection" : false, + "bfecc_number_substeps" : 10 }""") default_settings.AddMissingParameters(super(NavierStokesMPITwoFluidsSolver, cls).GetDefaultSettings()) @@ -81,156 +76,160 @@ def __init__(self, model, custom_settings): self.element_name = "TwoFluidNavierStokes" self.condition_name = "NavierStokesWallCondition" + self.element_integrates_in_time = True self.element_has_nodal_properties = True - self.min_buffer_size = 3 - if (self.settings["solver_type"].GetString() == "TwoFluids"): - self.element_name = "TwoFluidNavierStokes" + self.min_buffer_size = 3 - ## Construct the linear solver - self.trilinos_linear_solver = trilinos_linear_solver_factory.ConstructSolver(self.settings["linear_solver_settings"]) + self._bfecc_convection = self.settings["bfecc_convection"].GetBool() + if self._bfecc_convection: + self._bfecc_convection = False + KratosMultiphysics.Logger.PrintWarning(self.__class__.__name__, "BFECC is not implemented in MPI yet. Switching to standard level set convection.") - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPITwoFluidsSolver","Construction of NavierStokesMPITwoFluidsSolver finished.") + dynamic_tau = self.settings["formulation"]["dynamic_tau"].GetDouble() + self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, dynamic_tau) + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"Construction of NavierStokesMPITwoFluidsSolver finished.") def AddVariables(self): - super(NavierStokesMPITwoFluidsSolver, self).AddVariables() self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.PARTITION_INDEX) - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPITwoFluidsSolver","Variables for the Trilinos Two Fluid solver added correctly.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"Variables for the Trilinos Two Fluid solver added correctly.") def ImportModelPart(self): ## Construct the distributed import model part utility self.distributed_model_part_importer = DistributedImportModelPartUtility(self.main_model_part, self.settings) ## Execute the Metis partitioning and reading self.distributed_model_part_importer.ImportModelPart() - ## Sets DENSITY, VISCOSITY and SOUND_VELOCITY - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPITwoFluidsSolver","MPI model reading finished.") + KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__,"MPI model reading finished.") def PrepareModelPart(self): super(NavierStokesMPITwoFluidsSolver,self).PrepareModelPart() ## Construct the MPI communicators self.distributed_model_part_importer.CreateCommunicators() - def AddDofs(self): - super(NavierStokesMPITwoFluidsSolver, self).AddDofs() - - KratosMultiphysics.Logger.PrintInfo("NavierStokesMPITwoFluidsSolver","DOFs for the Trilinos Two Fluid solver added correctly.") - - - def Initialize(self): - - ## Construct the communicator - self.EpetraCommunicator = KratosTrilinos.CreateCommunicator() - - ## Get the computing model part - self.computing_model_part = self.GetComputingModelPart() - - KratosMultiphysics.NormalCalculationUtils().CalculateOnSimplex(self.computing_model_part, self.computing_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]) - - self.neighbour_search = KratosMultiphysics.FindNodalNeighboursProcess(self.computing_model_part) - (self.neighbour_search).Execute() - - self.accelerationLimitationUtility = KratosMultiphysics.FluidDynamicsApplication.AccelerationLimitationUtilities(self.computing_model_part, 5.0) - - ## If needed, create the estimate time step utility - if (self.settings["time_stepping"]["automatic_time_step"].GetBool()): - self.EstimateDeltaTimeUtility = self._GetAutomaticTimeSteppingUtility() - - # Set the time discretization utility to compute the BDF coefficients - time_order = self.settings["time_order"].GetInt() - if time_order == 2: - self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) + def _GetEpetraCommunicator(self): + if not hasattr(self, '_epetra_communicator'): + self._epetra_communicator = KratosTrilinos.CreateCommunicator() + return self._epetra_communicator + + def _CreateScheme(self): + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + # Cases in which the element manages the time integration + if self.element_integrates_in_time: + # "Fake" scheme for those cases in where the element manages the time integration + # It is required to perform the nodal update once the current time step is solved + scheme = KratosTrilinos.TrilinosResidualBasedIncrementalUpdateStaticSchemeSlip( + domain_size, + domain_size + 1) + # In case the BDF2 scheme is used inside the element, the BDF time discretization utility is required to update the BDF coefficients + if (self.settings["time_scheme"].GetString() == "bdf2"): + time_order = 2 + self.time_discretization = KratosMultiphysics.TimeDiscretization.BDF(time_order) + else: + err_msg = "Requested elemental time scheme \"" + self.settings["time_scheme"].GetString()+ "\" is not available.\n" + err_msg += "Available options are: \"bdf2\"" + raise Exception(err_msg) + # Cases in which a time scheme manages the time integration else: - raise Exception("Only \"time_order\" equal to 2 is supported. Provided \"time_order\": " + str(time_order)) - - ## Creating the Trilinos convergence criteria - self.conv_criteria = KratosTrilinos.TrilinosUPCriteria(self.settings["relative_velocity_tolerance"].GetDouble(), - self.settings["absolute_velocity_tolerance"].GetDouble(), - self.settings["relative_pressure_tolerance"].GetDouble(), - self.settings["absolute_pressure_tolerance"].GetDouble()) - - (self.conv_criteria).SetEchoLevel(self.settings["echo_level"].GetInt()) - - #### ADDING NEW PROCESSES : level-set-convection and variational-distance-process - self.level_set_convection_process = self._set_level_set_convection_process() - self.variational_distance_process = self._set_variational_distance_process() - - ## Creating the Trilinos incremental update time scheme (the time integration is defined within the TwoFluidNavierStokes element) - self.time_scheme = KratosTrilinos.TrilinosResidualBasedIncrementalUpdateStaticSchemeSlip(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE], # Domain size (2,3) - self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE]+1) # DOFs (3,4) - - - ## Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 3: + err_msg = "Custom scheme creation is not allowed. Two-fluids Navier-Stokes elements manage the time integration internally." + raise Exception(err_msg) + return scheme + + def _CreateLinearSolver(self): + linear_solver_configuration = self.settings["linear_solver_settings"] + return trilinos_linear_solver_factory.ConstructSolver(linear_solver_configuration) + + def _CreateConvergenceCriterion(self): + convergence_criterion = KratosTrilinos.TrilinosUPCriteria( + self.settings["relative_velocity_tolerance"].GetDouble(), + self.settings["absolute_velocity_tolerance"].GetDouble(), + self.settings["relative_pressure_tolerance"].GetDouble(), + self.settings["absolute_pressure_tolerance"].GetDouble()) + convergence_criterion.SetEchoLevel(self.settings["echo_level"].GetInt()) + return convergence_criterion + + def _CreateBuilderAndSolver(self): + # Set the guess_row_size (guess about the number of zero entries) for the Trilinos builder and solver + domain_size = self.GetComputingModelPart().ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + if domain_size == 3: guess_row_size = 20*4 - elif self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: + else: guess_row_size = 10*3 - - ## Construct the Trilinos builder and solver - if self.settings["consider_periodic_conditions"].GetBool() == True: - self.builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolverPeriodic(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver, - KratosFluid.PATCH_INDEX) + # Construct the Trilinos builder and solver + trilinos_linear_solver = self._GetLinearSolver() + epetra_communicator = self._GetEpetraCommunicator() + if self.settings["consider_periodic_conditions"].GetBool(): + builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolverPeriodic( + epetra_communicator, + guess_row_size, + trilinos_linear_solver, + KratosFluid.PATCH_INDEX) else: - self.builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolver(self.EpetraCommunicator, - guess_row_size, - self.trilinos_linear_solver) - - ## Construct the Trilinos Newton-Raphson strategy - self.solver = KratosTrilinos.TrilinosNewtonRaphsonStrategy(self.main_model_part, - self.time_scheme, - self.trilinos_linear_solver, - self.conv_criteria, - self.builder_and_solver, - self.settings["maximum_iterations"].GetInt(), - self.settings["compute_reactions"].GetBool(), - self.settings["reform_dofs_at_each_step"].GetBool(), - self.settings["move_mesh_flag"].GetBool()) - - (self.solver).SetEchoLevel(self.settings["echo_level"].GetInt()) - (self.solver).Initialize() - (self.solver).Check() - - self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.DYNAMIC_TAU, self.settings["formulation"]["dynamic_tau"].GetDouble()) - - - def _set_level_set_convection_process(self): + builder_and_solver = KratosTrilinos.TrilinosBlockBuilderAndSolver( + epetra_communicator, + guess_row_size, + trilinos_linear_solver) + return builder_and_solver + + def _CreateSolutionStrategy(self): + computing_model_part = self.GetComputingModelPart() + time_scheme = self._GetScheme() + linear_solver = self._GetLinearSolver() + convergence_criterion = self._GetConvergenceCriterion() + builder_and_solver = self._GetBuilderAndSolver() + return KratosTrilinos.TrilinosNewtonRaphsonStrategy( + computing_model_part, + time_scheme, + linear_solver, + convergence_criterion, + builder_and_solver, + self.settings["maximum_iterations"].GetInt(), + self.settings["compute_reactions"].GetBool(), + self.settings["reform_dofs_at_each_step"].GetBool(), + self.settings["move_mesh_flag"].GetBool()) + + def _CreateLevelSetConvectionProcess(self): # Construct the level set convection process - if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: + domain_size = self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] + linear_solver = self._GetLinearSolver() + computing_model_part = self.GetComputingModelPart() + epetra_communicator = self._GetEpetraCommunicator() + if domain_size == 2: level_set_convection_process = KratosTrilinos.TrilinosLevelSetConvectionProcess2D( - self.EpetraCommunicator, + epetra_communicator, KratosMultiphysics.DISTANCE, - self.computing_model_part, - self.trilinos_linear_solver) + computing_model_part, + linear_solver) else: level_set_convection_process = KratosTrilinos.TrilinosLevelSetConvectionProcess3D( - self.EpetraCommunicator, + epetra_communicator, KratosMultiphysics.DISTANCE, - self.computing_model_part, - self.trilinos_linear_solver) + computing_model_part, + linear_solver) return level_set_convection_process - - def _set_variational_distance_process(self): + def _CreateVariationalDistanceProcess(self): # Construct the variational distance calculation process - maximum_iterations = 2 + maximum_iterations = 2 #TODO: Make this user-definable + linear_solver = self._GetLinearSolver() + computing_model_part = self.GetComputingModelPart() + epetra_communicator = self._GetEpetraCommunicator() if self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2: variational_distance_process = KratosTrilinos.TrilinosVariationalDistanceCalculationProcess2D( - self.EpetraCommunicator, - self.computing_model_part, - self.trilinos_linear_solver, + epetra_communicator, + computing_model_part, + linear_solver, maximum_iterations, KratosMultiphysics.VariationalDistanceCalculationProcess2D.CALCULATE_EXACT_DISTANCES_TO_PLANE) else: variational_distance_process = KratosTrilinos.TrilinosVariationalDistanceCalculationProcess3D( - self.EpetraCommunicator, - self.computing_model_part, - self.trilinos_linear_solver, + epetra_communicator, + computing_model_part, + linear_solver, maximum_iterations, KratosMultiphysics.VariationalDistanceCalculationProcess3D.CALCULATE_EXACT_DISTANCES_TO_PLANE) diff --git a/applications/FluidDynamicsApplication/tests/AdjointVMSSensitivity2DTest/mpi_cylinder_test_parameters.json b/applications/FluidDynamicsApplication/tests/AdjointVMSSensitivity2DTest/mpi_cylinder_test_parameters.json index b67bdb0d5c70..a9577b304708 100644 --- a/applications/FluidDynamicsApplication/tests/AdjointVMSSensitivity2DTest/mpi_cylinder_test_parameters.json +++ b/applications/FluidDynamicsApplication/tests/AdjointVMSSensitivity2DTest/mpi_cylinder_test_parameters.json @@ -64,6 +64,11 @@ "element_type" : "vms", "use_orthogonal_subscales" : false, "dynamic_tau" : 1.0 + }, + "linear_solver_settings":{ + "solver_type" : "amgcl", + "max_iteration" : 200, + "tolerance" : 1e-8 } }, "processes": { diff --git a/applications/FluidDynamicsApplication/tests/darcy_channel_test.py b/applications/FluidDynamicsApplication/tests/darcy_channel_test.py index aa3b7fd2a0e0..578fcffc83eb 100644 --- a/applications/FluidDynamicsApplication/tests/darcy_channel_test.py +++ b/applications/FluidDynamicsApplication/tests/darcy_channel_test.py @@ -20,7 +20,7 @@ def InitializeSolutionStep(self): (self.time_discretization).ComputeAndSaveBDFCoefficients(self.GetComputingModelPart().ProcessInfo) # Initialize the solver current step - (self.solver).InitializeSolutionStep() + self._GetSolutionStrategy().InitializeSolutionStep() class DarcyChannelTest(UnitTest.TestCase): diff --git a/applications/StructuralMechanicsApplication/python_scripts/trilinos_structural_mechanics_solver.py b/applications/StructuralMechanicsApplication/python_scripts/trilinos_structural_mechanics_solver.py index 070b11f519a6..b90a7ab4d9b9 100644 --- a/applications/StructuralMechanicsApplication/python_scripts/trilinos_structural_mechanics_solver.py +++ b/applications/StructuralMechanicsApplication/python_scripts/trilinos_structural_mechanics_solver.py @@ -63,7 +63,7 @@ def Finalize(self): #### Specific internal functions #### - def get_epetra_communicator(self): + def _GetEpetraCommunicator(self): if not hasattr(self, '_epetra_communicator'): self._epetra_communicator = self._create_epetra_communicator() return self._epetra_communicator @@ -88,7 +88,7 @@ def _create_builder_and_solver(self): KratosMultiphysics.Logger.PrintWarning("Constraints are not yet implemented in MPI and will therefore not be considered!") linear_solver = self.get_linear_solver() - epetra_communicator = self.get_epetra_communicator() + epetra_communicator = self._GetEpetraCommunicator() if(self.main_model_part.ProcessInfo[KratosMultiphysics.DOMAIN_SIZE] == 2): guess_row_size = 15 else: