Nested Newton running for FSM, Arte 2D running, Arte 3D in preparation

ugbase/bridge/disc_bridges/algebra_bridge.cpp

@@ -66,6 +66,14 @@
 #include "lib_disc/operator/composite_conv_check.h"
 #include "lib_disc/spatial_disc/local_to_global/local_to_global_mapper.h"
 
+#include "lib_disc/operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h"
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+
+#include "lib_disc/operator/non_linear_operator/newton_solver/newtonUpdaterGeneric.h"
+
+#endif
+
 using namespace std;
 
 
@@ -265,6 +273,23 @@ static void Algebra(Registry& reg, string parentGroup)
 			.set_construct_as_smart_pointer(true);
 		reg.add_class_to_group(name, "AssembledLinearOperator", tag);
 	}
+
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+	//	generic Newton updater
+	{
+		std::string grp = parentGroup; grp.append("/Discretization");
+		using T = NewtonUpdaterGeneric<vector_type>;
+		string name = string("NewtonUpdaterGeneric").append(suffix);
+		reg.add_class_<T>(name, grp)
+			.add_constructor()
+//			.template add_constructor<void (*)(SmartPtr<NewtonUpdaterGeneric<vector_type> >)>("NewtonUpdaterGeneric")
+//			.add_method("updateNewton", static_cast<void(T::*)( vector_type &, vector_type const &, bool )> (&T::updateNewton))
+//			.add_method("updateNewton", static_cast<void(T::*)(vector_type &, number,  vector_type const &, number, vector_type const )> (&T::updateNewton))
+			.set_construct_as_smart_pointer(true);
+		reg.add_class_to_group(name, "NewtonUpdaterGeneric", tag);
+	}
+#endif
 
 
 //	NewtonSolver
@@ -301,6 +326,9 @@ static void Algebra(Registry& reg, string parentGroup)
 			.add_method("clear_inner_step_update", &T::clear_inner_step_update, "clear inner step update", "")
 			.add_method("add_step_update", &T::add_step_update, "data update called before every Newton step", "")
 			.add_method("clear_step_update", &T::clear_step_update, "clear step update", "")
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+			.add_method("setNewtonUpdater", &T::setNewtonUpdater, "set the Newton updater", "")
+#endif
 			.add_method("config_string", &T::config_string)
 			.set_construct_as_smart_pointer(true);
 		reg.add_class_to_group(name, "NewtonSolver", tag);

ugbase/bridge/grid_bridges/file_io_bridge.cpp

@@ -35,7 +35,6 @@
 #include "lib_grid/multi_grid.h"
 #include "lib_grid/file_io/file_io.h"
 #include "lib_grid/file_io/file_io_ugx.h"
-#include "lib_grid/file_io/file_io_vtu.h"
 
 using namespace std;
 
@@ -78,13 +77,6 @@ bool SaveGridHierarchy(MultiGrid& mg, const char* filename)
 	return SaveGridToFile(mg, mg.get_hierarchy_handler(), filename);
 }
 
-void SetVTURegionOfInterestIdentifier( char const * regOfInt )
-{
-	PROFILE_FUNC_GROUP("grid");
-	GridReaderVTU::setRegionOfInterestIdentifier( std::string( regOfInt ) );
-	return;
-}
-
 
 void RegisterGridBridge_FileIO(Registry& reg, string parentGroup)
 {
@@ -137,8 +129,7 @@ void RegisterGridBridge_FileIO(Registry& reg, string parentGroup)
 		.add_function("SaveParallelGridLayout", &SaveParallelGridLayout,
 				grp, "", "mg#filename#offset")
 		.add_function("SaveSurfaceViewTransformed", &SaveSurfaceViewTransformed)
-		.add_function("SaveGridLevelToFile", &SaveGridLevelToFile)
-		.add_function("SetVTURegionOfInterestIdentifier", static_cast<void (*)(char const *)>(&SetVTURegionOfInterestIdentifier) );
+		.add_function("SaveGridLevelToFile", &SaveGridLevelToFile);
 }
 
 }//	end of namespace

ugbase/lib_disc/CMakeLists.txt

@@ -34,6 +34,28 @@ project(P_LIB_DISCRETIZATION)
 
 include("../../cmake/ug_includes.cmake")
 
+option(ENABLE_NESTED_NEWTON "If enabled, nested Newton may be used, but needs for specification, else fallback to normal Newton" OFF)
+message(STATUS "")
+message(STATUS "Info: ${pluginName} options:")
+message(STATUS "    * ENABLE_NESTED_NEWTON:    ${ENABLE_NESTED_NEWTON} (options are: ON, OFF)")
+
+if( ENABLE_NESTED_NEWTON )
+
+set( SET_NESTED_NEWTON 1 ) 
+
+else(ENABLE_NESTED_NEWTON)
+
+set( SET_NESTED_NEWTON 0) 
+
+endif(ENABLE_NESTED_NEWTON)
+
+configure_file (
+  "${CMAKE_CURRENT_SOURCE_DIR}//operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h.in"
+  "${CMAKE_CURRENT_SOURCE_DIR}//operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h"
+  )
+
+
+
 set(srcDiscretization	domain.cpp
 						domain_util.cpp
 

ugbase/lib_disc/operator/non_linear_operator/line_search.h

@@ -38,9 +38,20 @@
 #include <vector>
 #include <cmath>
 #include <sstream>
+#include <limits>
 
 #include "common/common.h"
 
+#include "lib_disc/operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h"
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+
+#include "lib_disc/operator/non_linear_operator/newton_solver/newtonUpdaterGeneric.h"
+
+#endif
+
+
+
 namespace ug{
 
 ///////////////////////////////////////////////////////////
@@ -93,6 +104,10 @@ class ILineSearch
 
 	/// virtual destructor
 		virtual ~ILineSearch() {}
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+		virtual void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<TVector> > nU ) = 0;
+#endif
 };
 
 /// standard implementation of the line search based on the "sufficient descent"
@@ -108,18 +123,30 @@ class StandardLineSearch : public ILineSearch<TVector>
 		StandardLineSearch()
 		 :	 m_maxSteps(10), m_lambdaStart(1.0), m_lambdaReduce(0.5), m_alpha(0.25),
 		  	 m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(false), m_bCheckAll(false), m_offset("")
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+			 ,
+			m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
+#endif
 			 {};
 
 	///	constructor
 		StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest)
 		 :	 m_maxSteps(maxSteps), m_lambdaStart(lambdaStart), m_lambdaReduce(lambdaReduce), m_alpha(0.25),
 			 m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(bAcceptBest), m_bCheckAll(false), m_offset("")
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+			 ,
+			 m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
+#endif
 			 {};
 
 	///	constructor
 		StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest, bool bCheckAll)
 		 :	 m_maxSteps(maxSteps), m_lambdaStart(lambdaStart), m_lambdaReduce(lambdaReduce), m_alpha(0.25),
 			 m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(bAcceptBest), m_bCheckAll(bCheckAll), m_offset("")
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+			 ,
+			 m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
+#endif
 			 {};
 
 	///	returns information about configuration parameters
@@ -188,6 +215,37 @@ class StandardLineSearch : public ILineSearch<TVector>
 		//	loop line search steps
 			for(int k = 1; k <= m_maxSteps; ++k)
 			{
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+				// 	try on line u := u - lambda*p
+
+				//if( m_newtonUpdater != NULL )
+				//{
+				bool acceptedNewtonUpdate = m_newtonUpdater->updateSolution(u, 1.0, u, (-1)*lambda, p);
+
+//					if( ! m_newtonUpdater->updateNewton(u, 1.0, u, (-1)*lambda, p) )
+				if( ! acceptedNewtonUpdate )
+//					if( m_newtonUpdater->updateNewton(u,p) )
+				{
+					UG_LOG("Update in Line Search did not work.\n");
+					norm = std::numeric_limits<number>::max();
+						//VecScaleAdd(u, 1.0, u, (-1)*lambda, p);
+					//return false;
+				}
+				else
+				{
+					// 	compute new Defect
+					spOp->prepare(u);
+					spOp->apply(d, u);
+
+					//	compute new Residuum
+					norm = d.norm();
+				}
+				//}
+//				else
+//				{
+//					UG_LOG("Please set a Newton updater in Line search");
+//				}
+#else
 			// 	try on line u := u - lambda*p
 				VecScaleAdd(u, 1.0, u, (-1)*lambda, p);
 
@@ -198,6 +256,9 @@ class StandardLineSearch : public ILineSearch<TVector>
 			//	compute new Residuum
 				norm = d.norm();
 
+#endif
+
+
 			// 	compute reduction
 				vRho.push_back(norm/norm_old);
 
@@ -254,15 +315,42 @@ class StandardLineSearch : public ILineSearch<TVector>
 						UG_LOG(m_offset << "   ++++ Accept Best: Accepting step " <<
 							best+1 << ", Rate = "<< vRho[best] <<".\n");
 
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+					// 	try on line u := u - lambda*p
+
+					//if( m_newtonUpdater != NULL )
+					//{
+					if( ! m_newtonUpdater->updateSolution(u, 1.0, s, (-1)*m_lambdaStart*std::pow(m_lambdaReduce, (number)best), p) )
+					{
+						UG_LOG("Update in Line Search kmax did not work.\n");
+
+						norm = std::numeric_limits<number>::max();
+						//return false;
+					}
+					else
+					{
+						spOp->prepare(u);
+						spOp->apply(d, u);
+
+						// compute new Residuum
+						norm = d.norm();
+					}
+					//}
+//					else
+//					{
+//						UG_LOG("Please set a Newton updater in Line search");
+//					}
+#else
 				// 	try on line u := u - lambda*p
 					VecScaleAdd(u, 1.0, s, (-1)*m_lambdaStart*std::pow(m_lambdaReduce, (number)best), p);
-
 				// 	compute new Defect
 					spOp->prepare(u);
 					spOp->apply(d, u);
 
 					// compute new Residuum
 					norm = d.norm();
+#endif
+
 
 					// check if defect-norm is smaller than maximum allowed defect value
 					if (norm > m_maxDefect)
@@ -276,8 +364,15 @@ class StandardLineSearch : public ILineSearch<TVector>
 					break;
 				}
 
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+
+				// 	reset u and eventual local variables
+				m_newtonUpdater->resetSolution(u,s);
+#else
 			// 	reset u
 				u = s;
+#endif
+
 			}
 
 		//	print end line
@@ -288,10 +383,28 @@ class StandardLineSearch : public ILineSearch<TVector>
 					UG_LOG(m_offset << "   ++++ Line Search converged.\n");
 			}
 
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+			if( ! m_newtonUpdater->tellAndFixUpdateEvents(u) )
+			{
+				UG_LOG("unable to fix local Newton updates" << std::endl );
+				return false;
+			}
+#endif
 		//	we're done
 			return true;
 		}
 
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+
+		virtual void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<TVector> > nU )
+		{
+			m_newtonUpdater = nU;
+		}
+
+#endif
+
 	protected:
 	/// solution in line direction
 		vector_type s;
@@ -323,6 +436,11 @@ class StandardLineSearch : public ILineSearch<TVector>
 
 	/// number of spaces inserted before output
 		std::string m_offset;
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+private:
+		SmartPtr<NewtonUpdaterGeneric<TVector> > m_newtonUpdater;
+#endif
 };
 
 } // end namespace ug

ugbase/lib_disc/operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h.in

@@ -0,0 +1,17 @@
+/*
+ * nestedNewtonRFSwitch.h
+ *
+ *  Created on: 27.07.2021
+ *      Author: Markus Knodel
+ */
+
+#ifndef UGCORE_UGBASE_NESTEDNEWTONRFSWITCH_H_
+#define UGCORE_UGBASE_NESTEDNEWTONRFSWITCH_H_
+
+#ifndef ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+#define ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE ${SET_NESTED_NEWTON}
+#endif
+
+
+
+#endif /* UGCORE_UGBASE_NESTEDNEWTONRFSWITCH_H_ */

ugbase/lib_disc/operator/non_linear_operator/newton_solver/newton.h

@@ -47,6 +47,14 @@
 #include "newton_update_interface.h"
 #include "lib_algebra/operator/debug_writer.h"
 
+#include "nestedNewtonRFSwitch.h"
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+
+#include "newtonUpdaterGeneric.h"
+
+#endif
+
 namespace ug {
 
 /// Newton solver for assembling based discretizations
@@ -148,6 +156,13 @@ class NewtonSolver
 		void set_reassemble_J_freq(int freq)
 			{m_reassembe_J_freq = freq;};
 
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+		void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<vector_type> > nU )
+		{
+			m_newtonUpdater = nU;
+		}
+#endif
+
 	private:
 	///	help functions for debug output
 	///	\{
@@ -189,6 +204,13 @@ class NewtonSolver
 		std::vector<number> m_vNonLinSolverRates;
 		std::vector<number> m_vLinSolverRates;
 	/// \}
+
+#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
+
+		SmartPtr<NewtonUpdaterGeneric<vector_type> > m_newtonUpdater;
+
+#endif
+
 };
 
 }