diff --git a/classes/Integrator/EML_noCons.m b/classes/Integrator/EML_noCons.m index c954d50..831cf1b 100644 --- a/classes/Integrator/EML_noCons.m +++ b/classes/Integrator/EML_noCons.m @@ -35,7 +35,7 @@ end - function [resi, tang] = compute_resi_tang(self, zn1, zn, this_system) + function [resi, tang] = compute_resi_tang(self, zn1, zn, this_system, time_n) %% Abbreviations h = self.DT; diff --git a/classes/Integrator/EML_null.m b/classes/Integrator/EML_null.m index b756f0c..7c58e41 100644 --- a/classes/Integrator/EML_null.m +++ b/classes/Integrator/EML_null.m @@ -40,7 +40,7 @@ z0 = [this_simulation.Q_0',p0', this_simulation.V_0']; end - function [resi, tang] = compute_resi_tang(self, xn1, zn, this_system) + function [resi, tang] = compute_resi_tang(self, xn1, zn, this_system, time_n) %% Abbreviations h = self.DT; diff --git a/classes/Integrator/EML_reduced.m b/classes/Integrator/EML_reduced.m index 98bd7d4..33e2774 100644 --- a/classes/Integrator/EML_reduced.m +++ b/classes/Integrator/EML_reduced.m @@ -40,7 +40,7 @@ z0 = [this_simulation.Q_0',p0', this_simulation.V_0', self.LM0']; end - function [resi, tang] = compute_resi_tang(self, xn1, zn, this_system) + function [resi, tang] = compute_resi_tang(self, xn1, zn, this_system, time_n) %% Abbreviations h = self.DT; diff --git a/classes/Integrator/EMS_std.m b/classes/Integrator/EMS_std.m index 1338f5d..3bb0c98 100755 --- a/classes/Integrator/EMS_std.m +++ b/classes/Integrator/EMS_std.m @@ -40,7 +40,7 @@ end - function [resi, tang] = compute_resi_tang(self, zn1, zn, this_system) + function [resi, tang] = compute_resi_tang(self, zn1, zn, this_system, time_n) %% Abbreviations h = self.DT; @@ -55,24 +55,24 @@ lambda_n1 = zn1(2*nDOF+1:2*nDOF+m); Vext_n1 = this_system.external_potential(qn1); g_n1 = this_system.constraint(qn1); - %if ismethod(this_system,'get_inverse_mass_matrix') + if ismethod(this_system,'get_inverse_mass_matrix') IMn1 = this_system.get_inverse_mass_matrix(qn1); - %else - % Mn1 = this_system.get_mass_matrix(qn1); - % IMn1 = eye(size(Mn1)) / Mn1; - %end + else + Mn1 = this_system.get_mass_matrix(qn1); + IMn1 = eye(size(Mn1)) / Mn1; + end %% Known quantities from last time-step qn = zn(1:nDOF); pn = zn(nDOF+1:2*nDOF); Vext_n = this_system.external_potential(qn); - %if ismethod(this_system,'get_inverse_mass_matrix') + if ismethod(this_system,'get_inverse_mass_matrix') IMn = this_system.get_inverse_mass_matrix(qn); - %else - % Mn = this_system.get_mass_matrix(qn); - % IMn = eye(size(Mn)) / Mn; - %end + else + Mn = this_system.get_mass_matrix(qn); + IMn = eye(size(Mn)) / Mn; + end %% MP evaluated quantities q_n05 = 0.5 * (qn + qn1); @@ -80,12 +80,12 @@ DVext_n05 = this_system.external_potential_gradient(q_n05); D_1_T_n05 = this_system.kinetic_energy_gradient_from_momentum(q_n05, p_n05); - %if ismethod(this_system,'get_inverse_mass_matrix') + if ismethod(this_system,'get_inverse_mass_matrix') IMn05 = this_system.get_inverse_mass_matrix(q_n05); - %else - % Mn05 = this_system.get_mass_matrix(q_n05); - % IMn05 = eye(size(Mn05)) / Mn05; - %end + else + Mn05 = this_system.get_mass_matrix(q_n05); + IMn05 = eye(size(Mn05)) / Mn05; + end % kinetic energy with mixed evaluations T_qn1pn = 0.5 * pn' * IMn1 * pn; diff --git a/classes/Integrator/MP_std.m b/classes/Integrator/MP_std.m index 496f90f..8f56cc1 100755 --- a/classes/Integrator/MP_std.m +++ b/classes/Integrator/MP_std.m @@ -34,7 +34,7 @@ end - function [resi, tang] = compute_resi_tang(self, zn1, zn, this_system) + function [resi, tang] = compute_resi_tang(self, zn1, zn, this_system, time_n) %% Abbreviations M = this_system.MASS_MAT; diff --git a/classes/Metis/Metis.m b/classes/Metis/Metis.m index 5e5570e..08c4f4c 100755 --- a/classes/Metis/Metis.m +++ b/classes/Metis/Metis.m @@ -27,6 +27,7 @@ NUM_ITER MEAN_ITER ALPHA_0 + ADVANCED_INIT = false; %% Postprocessing parameters % Can be given in an input-file (not necessary for computing) diff --git a/classes/Solver/Solver.m b/classes/Solver/Solver.m index 59362aa..5867eca 100755 --- a/classes/Solver/Solver.m +++ b/classes/Solver/Solver.m @@ -5,6 +5,7 @@ properties MAX_ITERATIONS TOLERANCE + ADVANCED_INIT end methods @@ -15,6 +16,7 @@ % Set Newton tolerance and max. iterations from simulation self.TOLERANCE = this_simulation.TOLERANCE; self.MAX_ITERATIONS = this_simulation.MAX_ITERATIONS; + self.ADVANCED_INIT = this_simulation.ADVANCED_INIT; end @@ -50,7 +52,9 @@ zn = zn1; % Initial Guess - %zn1 = self.newton_initial_guess(zn1,this_system,this_simulation, this_integrator.DT); + if self.ADVANCED_INIT + zn1 = self.newton_initial_guess(zn1,this_system,this_simulation, this_integrator.DT); + end % Print current time-step fprintf(this_simulation.log_file_ID, '%s: %s\n', datestr(now, 0),'----------------------------------------------------'); diff --git a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/EMH_ref_sol.mat b/input/published/MUBO_kinon_betsch_2024/quat_freeRB/EMH_ref_sol.mat deleted file mode 100644 index 0f26377..0000000 Binary files a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/EMH_ref_sol.mat and /dev/null differ diff --git a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/EML_ref_sol.mat b/input/published/MUBO_kinon_betsch_2024/quat_freeRB/EML_ref_sol.mat deleted file mode 100644 index ba80c59..0000000 Binary files a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/EML_ref_sol.mat and /dev/null differ diff --git a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/error_analysis_rigidBodyRotatingQuat.m b/input/published/MUBO_kinon_betsch_2024/quat_freeRB/error_analysis_rigidBodyRotatingQuat.m deleted file mode 100755 index 6dd63a3..0000000 --- a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/error_analysis_rigidBodyRotatingQuat.m +++ /dev/null @@ -1,71 +0,0 @@ -%% System Parameters -% Name of system in /classes/System -%SYSTEM = 'RigidBodyRotatingQuaternions'; -SYSTEM = 'RigidBodyRotatingQuaternionsRegularMassMatrix'; -% External acceleration -EXT_ACC = 0; -% Initial configuration -Q_0 = [1; 0; 0; 0]; -% Initial velocity -Omega_0 = [10;20;20]; -%extract vector and scalar part form quaternion -Q0_vec = Q_0(2:4); -Q0_scalar = Q_0(1); - -%skew-sym matrix corresponding to vector part -Q0_hat = [0, -Q0_vec(3), Q0_vec(2); - Q0_vec(3), 0, -Q0_vec(1); - -Q0_vec(2), Q0_vec(1), 0]; - -% transformation matrix -G_Q0 = [-Q0_vec, Q0_scalar*eye(3) - Q0_hat]; -V_0 = 1/2*G_Q0'*Omega_0; -% Mass -MASS = 1; -% Spatial dimensions -DIM = 3; - -% clear unnecessary variables (crucial for further processing!) -clear Omega_0 Q0_hat Q0_vec Q0_scalar G_Q0 - -%% Integrator -% Name of routines in /classes/Integrator to be analyzed -%INTEGRATOR = {'EML'}; -INTEGRATOR = {'EMS_std'}; -%INTEGRATOR = {'MP_Livens'}; -%INTEGRATOR = {'MP_std'}; -% Parameters of the method -INT_PARA = [NaN, NaN]; -% time step sizes to be analyzed -DT = [0.01, 0.005, 0.001, 0.0005, 0.0001, 0.00001]; -% quantity which is to be analyzed -CONV_QUANTITY = 'q'; -% starting time -T_0 = 0; -% end time -T_END = 0.1; - -%% Solver Method -% maximum number of iterations of Newton Rhapson method -MAX_ITERATIONS = 40; -% tolerance of Newton Rhapson method -TOLERANCE = 1E-09; - -%% Postprocessing -% Animation of trajectory [true/false] -shouldAnimate = false; -% List of desired quantities for plotting in postprocessing -plot_quantities = {}; -% Export of simulation results in a .mat-file [true/false] -should_export = true; -% Export of figures in .eps- and .tikz-files -should_export_figures = true; -% Path where export-folder is created -export_path = 'scratch/'; -% Matlab2Tikz (metis searches for matlab2tikz here. if not available, it -% clones the matlab2tikz repository there) -matlab2tikz_directory = '~/git/matlab2tikz'; - -%% Write variables into a .mat-File -% for further processing by metis -save(mfilename); \ No newline at end of file diff --git a/input/published/MUBO_kinon_betsch_2024/quat_heavytop/error_analysis_HeavyTopQuat.m b/input/published/MUBO_kinon_betsch_2024/quat_heavytop/error_analysis_HeavyTopQuat.m index 77dafc4..e2e12af 100644 --- a/input/published/MUBO_kinon_betsch_2024/quat_heavytop/error_analysis_HeavyTopQuat.m +++ b/input/published/MUBO_kinon_betsch_2024/quat_heavytop/error_analysis_HeavyTopQuat.m @@ -1,7 +1,8 @@ %% System Parameters % Name of system in /classes/System -%SYSTEM = 'HeavyTopQuaternions'; -SYSTEM = 'HeavyTopQuaternionsRegularMassMatrix'; +SYSTEM = 'HeavyTopQuaternions'; %for EML,MPL +%SYSTEM = 'HeavyTopQuaternionsRegularMassMatrix'; % for EMH,MPH + % External acceleration g = 9.81; EXT_ACC = [0; 0; -g]; @@ -46,11 +47,10 @@ %% Integrator % Name of routines in /classes/Integrator to be analyzed -%INTEGRATOR = {'EML'}; %EMS_std -%INTEGRATOR = {'MP_Livens'}; -INTEGRATOR = {'MP_std'}; +INTEGRATOR = {'EML','MP_Livens'}; +%INTEGRATOR = {'EMS_std','MP_std'}; % Parameters of the method -INT_PARA = [NaN, NaN]; +INT_PARA = [NaN, NaN; NaN, NaN]; % time step sizes to be analyzed DT = [0.01, 0.005, 0.001, 0.0005]; % quantity which is to be analyzed @@ -65,6 +65,8 @@ MAX_ITERATIONS = 40; % tolerance of Newton Rhapson method TOLERANCE = 1E-09; +% advanced initial guess +ADVANCED_INIT = true; %true required for EMH and MPH %% Postprocessing % Animation of trajectory [true/false] diff --git a/input/published/MUBO_kinon_betsch_2024/quat_heavytop/single_analysis_HeavyTopQuat.m b/input/published/MUBO_kinon_betsch_2024/quat_heavytop/single_analysis_HeavyTopQuat.m index 65a660d..4db1122 100644 --- a/input/published/MUBO_kinon_betsch_2024/quat_heavytop/single_analysis_HeavyTopQuat.m +++ b/input/published/MUBO_kinon_betsch_2024/quat_heavytop/single_analysis_HeavyTopQuat.m @@ -1,7 +1,7 @@ %% System Parameters % Name of system in /classes/System SYSTEM = 'HeavyTopQuaternions'; -%SYSTEM = 'HeavyTopQuaternionsRegularMassMatrix'; +%SYSTEM = 'HeavyTopQuaternionsRegularMassMatrix'; %required for EMH, MPH % External acceleration g = 9.81; @@ -48,9 +48,9 @@ %% Integrator % Name of routine in /classes/Integrator INTEGRATOR = 'EML'; -%INTEGRATOR = 'EMS_std'; +%INTEGRATOR = 'EMS_std'; %is EMH %INTEGRATOR = 'MP_Livens'; -%INTEGRATOR = 'MP_std'; +%INTEGRATOR = 'MP_std'; %is MPH % Parameters of the method INT_PARA = [NaN, NaN]; % time step size @@ -65,6 +65,8 @@ MAX_ITERATIONS = 40; % tolerance of Newton Rhapson method TOLERANCE = 1E-09; +% advanced initial guess +ADVANCED_INIT = false; %true required for EMH and MPH %% Postprocessing % Animation of trajectory [true/false] diff --git a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/single_analysis_rigidBodyRotatingQuat.m b/input/published/MUBO_kinon_betsch_2024/single_analysis_rigidBodyRotatingQuat.m similarity index 98% rename from input/published/MUBO_kinon_betsch_2024/quat_freeRB/single_analysis_rigidBodyRotatingQuat.m rename to input/published/MUBO_kinon_betsch_2024/single_analysis_rigidBodyRotatingQuat.m index 536e0df..f2c3ce3 100644 --- a/input/published/MUBO_kinon_betsch_2024/quat_freeRB/single_analysis_rigidBodyRotatingQuat.m +++ b/input/published/MUBO_kinon_betsch_2024/single_analysis_rigidBodyRotatingQuat.m @@ -39,11 +39,11 @@ % Parameters of the method INT_PARA = [NaN, NaN]; % time step size -DT = 0.01; +DT = 0.05; % starting time T_0 = 0; % end time -T_END = 10; %2 +T_END = 2; %10; %% Solver Method % maximum number of iterations of Newton Rhapson method diff --git a/start_metis_error_analysis.m b/start_metis_error_analysis.m index 836cd43..3b968e0 100755 --- a/start_metis_error_analysis.m +++ b/start_metis_error_analysis.m @@ -24,7 +24,6 @@ % Metis creates a dummy simulation object, the system and solver from input-file [dummy_simulation, system, dummy_integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/quat_heavytop/error_analysis_HeavyTopQuat', 1, 1); -%[dummy_simulation, system, dummy_integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/quat_freeRB/error_analysis_rigidBodyRotatingQuat', 1, 1); % Check how many different timestepsizes and integrators are analyzed n_DT = numel(dummy_simulation.ALL_DT); n_INT = numel(dummy_simulation.ALL_INTEGRATOR); @@ -39,7 +38,6 @@ % Metis creates objects for current timestepsize and integrator [current_simulation, ~, current_integrator, ~] = Metis('input/published/MUBO_kinon_betsch_2024/quat_heavytop/error_analysis_HeavyTopQuat', i, j); - %[current_simulation, ~, current_integrator, ~] = Metis('input/published/MUBO_kinon_betsch_2024/quat_freeRB/error_analysis_rigidBodyRotatingQuat', i, j); %% METIS solver % Solve system with solver and current integrator diff --git a/start_metis_single_analysis.m b/start_metis_single_analysis.m index 3af28fb..7fa9e93 100755 --- a/start_metis_single_analysis.m +++ b/start_metis_single_analysis.m @@ -25,10 +25,11 @@ % Add all subdirectories and to the current path addpath(genpath(fileparts(which(mfilename)))); % Metis creates objects from input-file -%[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/single_analysis_RedundantMassSpring', 1, 1); +[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/single_analysis_RedundantMassSpring', 1, 1); %[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/single_analysis_SpringPendulumPolar', 1, 1); -%[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/quat_freeRB/single_analysis_rigidBodyRotatingQuat', 1, 1); -[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/single_analysis_ClosedLoopMBS', 1, 1); +%[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/single_analysis_rigidBodyRotatingQuat', 1, 1); +%[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/quat_heavytop/single_analysis_HeavyTopQuat', 1, 1); +%[simulation, system, integrator, solver] = Metis('input/published/MUBO_kinon_betsch_2024/single_analysis_ClosedLoopMBS', 1, 1); %% METIS solver % Solve system with chosen solver and integration scheme