last minor changes for publication

plkinon · Jun 17, 2024 · f3237b1 · f3237b1
1 parent 7019e3c
commit f3237b1
Show file tree

Hide file tree

Showing 15 changed files with 45 additions and 108 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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,37 +55,37 @@
             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);
             p_n05 = 0.5 * (pn + pn1);
             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
@@ -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
@@ -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
@@ -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
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
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
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
@@ -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
@@ -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/...B/single_analysis_rigidBodyRotatingQuat.m → ...4/single_analysis_rigidBodyRotatingQuat.m b/...B/single_analysis_rigidBodyRotatingQuat.m → ...4/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
@@ -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
@@ -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