Newton finally working with extreme underrelaxation

mfds/active_set_operator.cpp

@@ -165,7 +165,7 @@ namespace mfds {
     }
 
     Eigen::VectorXd heat_analytic_linear_operator_t::mult_internal(const Eigen::VectorXd &v) const {
-        std::cout << "heat: mult internal " << std::endl << std::flush;
+        //std::cout << "heat: mult internal " << std::endl << std::flush;
         auto const vb = vtint::split_to_block_vector(v, t_points.size() - 1);
         auto const res = this->operator*(vb);
         return res.to_vector();

mfds/active_set_operator.hpp

@@ -17,10 +17,10 @@
 #include <optional>
 
 
-
 namespace mfds {
     template<class Op_t>
     class boundary_condition_wrapper_t;
+
     class heat_analytic_linear_operator_t;
 }
 using Eigen::SparseMatrix;
@@ -68,7 +68,8 @@ namespace mfds {
         }
 
     public:
-        boundary_condition_wrapper_t(Op_t& Op_, std::vector<int> const& boundary_dofs_) : Op(Op_), boundary_dofs(boundary_dofs_){
+        boundary_condition_wrapper_t(Op_t &Op_, std::vector<int> const &boundary_dofs_) : Op(Op_), boundary_dofs(
+                boundary_dofs_) {
             is_boundary = std::vector<bool>(Op.cols(), false);
             for (auto const dof: boundary_dofs) {
                 is_boundary[dof] = true;
@@ -92,7 +93,7 @@ namespace mfds {
         }
 
         [[nodiscard]] Eigen::VectorXd mult_internal(Eigen::VectorXd const &v) const {
-            std::cout << "Bwrap mult" << std::endl << std::flush;
+            //std::cout << "Bwrap mult" << std::endl << std::flush;
             Eigen::VectorXd res = Op * v;
             for (auto const &dof: boundary_dofs) {
                 res(dof) = v(dof);
@@ -163,6 +164,148 @@ namespace mfds {
 
     };
 
+    template<class Op_t>
+    class active_set_wrapper_t{
+    public:
+        active_set_wrapper_t(Op_t &Op_, Eigen::VectorXd const &lower_bound_, Eigen::VectorXd const &upper_bound_,
+                             double const c_ = 1.)
+                : Op(Op_), lower_bound(lower_bound_), upper_bound(upper_bound_), c(c_) {
+            if (lower_bound.size() != Op.cols()) {
+                throw std::runtime_error("Lower bound size does not match operator size");
+            }
+            if (upper_bound.size() != Op.cols()) {
+                throw std::runtime_error("Upper bound size does not match operator size");
+            }
+            is_active.resize(Op.cols());
+            is_active.setConstant(false);
+            is_active_lower = is_active;
+            is_active_upper = is_active;
+
+        };
+
+        std::tuple<Eigen::VectorXd, Eigen::VectorXd>
+        solve(Eigen::VectorXd const &rhs, Eigen::VectorXd const &u0, Eigen::VectorXd const &lambda0,
+              double const omega, int const maxiter = 100) {
+
+            Eigen::VectorXd u = u0;
+            Eigen::VectorXd lambdas = lambda0;
+
+
+            Eigen::Index n_changed = u.size();
+            int iter = 0;
+            double u_diff = 1e10;
+            while ((n_changed > 0 || u_diff > 1e-3)&& iter < maxiter) {
+                auto const u_old = u;
+                auto const lambdas_old = lambdas;
+                auto const is_active_old = is_active;
+
+                //determine active sets
+                std::tie(u, lambdas) = determine_acitve_set(u, lambdas);
+
+                //now we can use the boundary condition wrapper to set ones on the diagonals of active sets
+                boundary_condition_wrapper_t K(Op, is_active_index);
+                Eigen::VectorXd cur_rhs = adapt_rhs_for_active_set(rhs, u);
+
+                //now we can solve the equation system using Eigen CG
+                Eigen::ConjugateGradient<decltype(K),
+                        Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> cg;
+                cg.compute(K);
+                Eigen::VectorXd const sol = cg.solve(cur_rhs);
+                std::tie(u, lambdas) = distribute_values(sol, u, lambdas);
+                n_changed = 0;
+                for (Eigen::Index i = 0; i < u.size(); ++i) {
+                    if (is_active_old[i] != is_active[i]) {
+                        n_changed++;
+                    }
+                }
+
+                if(iter == 0){
+                    n_changed = u.size();
+                }
+                iter++;
+                std::cout << "iter: " << iter << " n_changed: " << n_changed << std::endl;
+
+                u = omega*u + (1-omega)*u_old;
+                lambdas = omega*lambdas + (1-omega)*lambdas_old;
+
+                u_diff = (u-u_old).lpNorm<Eigen::Infinity>();
+                std::cout << "change in ||u-u_old||_infty: " << u_diff << std::endl;
+            }
+
+
+            return std::make_tuple(u, lambdas);
+
+        }
+
+    private:
+        std::tuple<Eigen::VectorXd, Eigen::VectorXd>
+        determine_acitve_set(Eigen::VectorXd const &u, Eigen::VectorXd const &lambdas) {
+            Eigen::VectorXd u_next = u;
+            Eigen::VectorXd lambda_next = lambdas;
+
+            is_active.setConstant(false);
+            is_active_lower.setConstant(false);
+            is_active_upper.setConstant(false);
+            is_active_index.clear();
+            for (Eigen::Index i = 0; i < u.size(); ++i) {
+                double const ui = u[i];
+                double const li = lambdas[i];
+                double const up = upper_bound[i];
+                double const um = lower_bound[i];
+                if (li + c * (um - ui) > 0) {
+                    u_next[i] = um;
+                    is_active_lower[i] = true;
+                } else if (li + c * (up - ui) < 0) {
+                    u_next[i] = up;
+                    is_active_upper[i] = true;
+                } else {
+                    lambda_next[i] = 0;
+                }
+                is_active[i] = is_active_lower[i] || is_active_upper[i];
+                if (is_active[i]) {
+                    is_active_index.push_back(i);
+                }
+            }
+            return std::make_tuple(u_next, lambda_next);
+        }
+
+        Eigen::VectorXd adapt_rhs_for_active_set(Eigen::VectorXd const &rhs, Eigen::VectorXd const &u) {
+            Eigen::VectorXd u_adapted = u;
+            for (Eigen::Index i = 0; i < u.size(); ++i) {
+                if (!is_active[i]) {
+                    u_adapted[i] = 0.;
+                }
+            }
+            Eigen::VectorXd const rhs_adapted = rhs - Op * u_adapted;
+            //std::cout << "rhs_adapted to rhs diff" << (rhs_adapted-rhs).lpNorm<Eigen::Infinity>() << std::endl;
+            return rhs_adapted;
+        }
+
+        std::tuple<Eigen::VectorXd, Eigen::VectorXd>
+        distribute_values(Eigen::VectorXd const &sol, Eigen::VectorXd const &u, Eigen::VectorXd const &lambda) {
+            Eigen::VectorXd u_next = u;
+            Eigen::VectorXd lambda_next = lambda;
+            for (Eigen::Index i = 0; i < u.size(); ++i) {
+                if (is_active[i]) {
+                    lambda_next[i] = -sol[i];
+                } else {
+                    u_next[i] = sol[i];
+                }
+            }
+            return std::make_tuple(u_next, lambda_next);
+        }
+
+    private:
+        Op_t &Op;
+        Eigen::VectorXd const lower_bound;
+        Eigen::VectorXd const upper_bound;
+        Eigen::ArrayX<bool> is_active;
+        Eigen::ArrayX<bool> is_active_lower;
+        Eigen::ArrayX<bool> is_active_upper;
+        std::vector<int> is_active_index;
+        double c;
+    };
+
 } // mfds
 
 
@@ -172,6 +315,7 @@ namespace mfds {
 namespace Eigen::internal {
 
     using namespace mfds;
+
     template<typename Rhs>
     struct generic_product_impl<heat_analytic_linear_operator_t, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
             : generic_product_impl_base<heat_analytic_linear_operator_t, Rhs, generic_product_impl<heat_analytic_linear_operator_t, Rhs> > {

playground/CMakeLists.txt

@@ -48,6 +48,12 @@ target_include_directories(example_matrix_free_operator_2d
         PUBLIC ${MFEM_INCLUDE_DIRS})
 target_compile_features(example_matrix_free_operator_2d PRIVATE cxx_std_17)
 
+add_executable(example_active_set_2d example_active_set_2d.cpp)
+TARGET_LINK_LIBRARIES(example_active_set_2d PUBLIC ${MFEM_LIBRARIES} hilbert_transformation_lib mfds vtint)
+target_include_directories(example_active_set_2d
+        PUBLIC ${MFEM_INCLUDE_DIRS})
+target_compile_features(example_matrix_free_operator_2d PRIVATE cxx_std_17)
+
 add_executable(mixed_iga_mfem mixed_iga_mfem.cpp)
 TARGET_LINK_LIBRARIES(mixed_iga_mfem PUBLIC ${MFEM_LIBRARIES} mfds vtint)
 target_include_directories(mixed_iga_mfem

playground/example_active_set_2d.cpp

+//
+// Created by mreic on 15.07.2024.
+//
+
+
+#include <iostream>
+
+#include "mfem.hpp"
+#include "active_set_operator.hpp"
+#include "heat_control.h"
+#include "target.hpp"
+#include "conjugate_gradient.hpp"
+#include "homogenization.hpp"
+#include "time_dependent_gridfunction.hpp"
+#include <Eigen/IterativeLinearSolvers>
+
+static double target_function(mfem::Vector const &p, double const t) {
+    double const x = p(0);
+    double const y = p(1);
+    //return 16*t * t * x * (1 - x) * y * (1 - y);
+    double const offset = 0.;
+    if(t<offset){
+        return 0.;
+    }
+    return std::sin(M_PI * x) * std::sin(M_PI * y) * std::sin(M_PI * (t-offset));
+}
+
+static double lower_bound_fun(double const t, mfem::Vector const &p) {
+    return 0.;
+}
+
+static double upper_bound_fun(double const t, mfem::Vector const &p) {
+    return .9;//std::min(target_function(p, t), .5);
+}
+
+int main(int argc, char *argv[]) {
+    int const refinement = 4;
+    int const nx = 2 * std::pow(2, refinement);
+    int const ny = 2 * std::pow(2, refinement);
+    int const nt = 2 * std::pow(2, refinement);
+    double const T = 1.;
+    double const cg_tol = 1e-12;
+    int const cg_max_iter = 500;
+    double const newton_c = 1;
+    double const newton_omega = 0.1;
+
+
+    mfem::Mesh spatial_mesh = mfem::Mesh::MakeCartesian2D(nx, ny, mfem::Element::TRIANGLE, true, 1.0, 1.0);
+    int const dimX = spatial_mesh.Dimension();
+    mfem::Mesh temporal_mesh = mfem::Mesh::MakeCartesian1D(nt, T);
+
+    //construct FE spaces
+    mfem::H1_FECollection X_fec(1, spatial_mesh.Dimension());
+    mfem::FiniteElementSpace X(&spatial_mesh, &X_fec);
+    mfem::H1_FECollection T_fec(1, temporal_mesh.Dimension());
+    mfem::FiniteElementSpace T_test(&temporal_mesh, &T_fec);
+
+    //construct matrix free operator
+    double const rho = std::pow(std::min(1. / nx, 1. / ny), 2);
+    mfds::heat_analytic_linear_operator_t heat_operator(X, nt, T, rho);
+
+
+    //get space time boundary
+    auto const space_time_boundary_dofs = heat_operator.get_space_time_boundary_dofs();
+    std::cout << "n_stb: " << space_time_boundary_dofs.size() << std::endl;
+    std::vector<bool> isBoundary(X.GetTrueVSize() * nt, false);
+    std::cout << "Space-Time Boundary dofs:" << std::endl;
+    for (auto const dof: space_time_boundary_dofs) {
+        std::cout << dof << std::endl;
+        isBoundary[dof] = true;
+    }
+
+    //calculate inhomogeneous rhs for target
+    auto const rhs_inhom = mfds::compute_spatially_inhom_rhs(target_function, T_test, X);
+
+    Eigen::VectorXd const u_bdr = Eigen::VectorXd::Zero(X.GetTrueVSize() * nt);
+
+    //create bdr condition wrapper
+    mfds::boundary_condition_wrapper_t boundary_condition_wrapper(heat_operator, space_time_boundary_dofs);
+    auto const rhs_bdr = boundary_condition_wrapper.modify_rhs(rhs_inhom.to_vector(), u_bdr);
+
+
+
+    //project lower and upper bound to gridfunction
+    vtint::time_dependent_gridfunction_t lower_bound_gf(T_test, X, lower_bound_fun);
+    vtint::time_dependent_gridfunction_t upper_bound_gf(T_test, X, upper_bound_fun);
+    Eigen::VectorXd lower_bound = vtint::remove_initial_condition(lower_bound_gf.get_blockvector()).to_vector();
+    Eigen::VectorXd upper_bound = vtint::remove_initial_condition(upper_bound_gf.get_blockvector()).to_vector();
+
+
+    //apply active set solver to solve problem
+    mfds::active_set_wrapper_t active_set_wrapper(boundary_condition_wrapper, lower_bound, upper_bound, newton_c);
+
+    //determine initial values
+    Eigen::VectorXd u0 = 0.5*(lower_bound+upper_bound);
+    Eigen::VectorXd lambda0 = -rhs_bdr + boundary_condition_wrapper*u0;
+
+    for(auto const& i : space_time_boundary_dofs){
+        u0(i) = 0.;
+        lambda0(i) = 0.;
+    }
+
+    auto [u_vec, lambda_vec] = active_set_wrapper.solve(rhs_bdr, u0, lambda0, newton_omega);
+
+
+    //add initial condition
+    auto const state_block = vtint::split_to_block_vector(u_vec, nt);
+    mfem::Vector u0_vec(X.GetTrueVSize());
+    u0_vec = 0.;
+    auto const state_full = vtint::add_initial_condition(state_block, u0_vec);
+
+    //compute L2 error, to be sure everything is correct
+    vtint::time_dependent_gridfunction_t state_gf(T_test, X, state_full);
+    //state_gf.save_as_vtk("state", 0.01);
+    auto target_fun_reverse_arguments = [](double const t, mfem::Vector const &p) {
+        return target_function(p, t);
+    };
+    double const target_state_l2_error = state_gf.l2_error(target_fun_reverse_arguments);
+
+    std::cout << "L2 error: " << target_state_l2_error << std::endl;
+
+    //save function to paraview
+    state_gf.save_as_vtk("state", 0.01);
+
+
+    vtint::time_dependent_gridfunction_t target_gf(T_test, X, target_fun_reverse_arguments);
+    target_gf.save_as_vtk("target", 0.01);
+
+    return 0;
+}
\ No newline at end of file

playground/example_matrix_free_operator_2d.cpp

@@ -55,7 +55,7 @@ int main(int argc, char *argv[]) {
     mfem::FiniteElementSpace T_test(&temporal_mesh, &T_fec);
 
     //construct matrix free operator
-    double const rho = 0*std::pow(std::min(1. / nx, 1. / ny), 2);
+    double const rho = std::pow(std::min(1. / nx, 1. / ny), 2);
     mfds::heat_analytic_linear_operator_t heat_operator(X, nt, T, rho);
 
 
@@ -91,19 +91,17 @@ int main(int argc, char *argv[]) {
     auto const [sim_res, sim_timing] = mfds::solve_heat_analytic(X, T_test, target_function, rho, true);
     auto const state_reference = sim_res.state_full;
 
-    /*
+
 
     //wrap operator so that boundary conditions can be applied
-    //with zero boundary condition
-    std::vector<double> const bVals(X.GetTrueVSize() * nt, 0.);
-    mfds::boundary_condition_wrapper_t boundary_condition_wrapper(&heat_operator, isBoundary, bVals);
-    auto const rhs_boundary_adapted = boundary_condition_wrapper.adapt_rhs(rhs_inhom.to_vector());
+
 
     //project lower and upper bound to gridfunction
     vtint::time_dependent_gridfunction_t lower_bound_gf(T_test, X, lower_bound_fun);
     vtint::time_dependent_gridfunction_t upper_bound_gf(T_test, X, upper_bound_fun);
     Eigen::VectorXd lower_bound = vtint::remove_initial_condition(lower_bound_gf.get_blockvector()).to_vector();
     Eigen::VectorXd upper_bound = vtint::remove_initial_condition(upper_bound_gf.get_blockvector()).to_vector();
+    /*
     mfds::active_set_solver_t active_set_solver(&boundary_condition_wrapper, rhs_boundary_adapted, isBoundary,
                                                 lower_bound,
                                                 upper_bound);