Skip to content
Snippets Groups Projects
Commit fbf503e8 authored by Reichelt, Michael's avatar Reichelt, Michael
Browse files

Still not converging...

parent 46ccb9ee
No related branches found
No related tags found
No related merge requests found
Pipeline #449155 failed
Show whitespace changes
Inline Side-by-side
mfds/active_set_operator.cpp
+ 16
− 0
View file @ fbf503e8
......@@ -52,6 +52,18 @@ namespace mfds {
generalized_eigenvalues = std::get<0>(values_and_vectors);
//get space time boundary dofs
get_space_time_boundary_dofs_internal(X, nt, T);
//get diagonal of mass
mass_diagonal = Eigen::VectorXd::Zero(X.GetVSize()*nt);
Eigen::Index cur_ind = 0;
for(int i=0; i<nt; ++i){
double const temp_val = mass_t.coeff(i, i);
for(int j=0; j<X.GetVSize(); ++j){
double const spatial_val = mass_x.coeff(j, j);
mass_diagonal[cur_ind] = temp_val*spatial_val;
cur_ind++;
}
}
}
void
......@@ -114,4 +126,8 @@ namespace mfds {
return 1.;
}
Eigen::VectorXd const &heat_analytic_linear_operator_t::get_diagonal_preconditioner() const {
return mass_diagonal;
}
} // mfds
\ No newline at end of file
mfds/active_set_operator.hpp
+ 27
− 7
View file @ fbf503e8
......@@ -68,14 +68,14 @@ namespace mfds {
Eigen::VectorXd solve(double const cg_tol, int const cg_max_iter, double const alpha = 0.5,
std::optional<Eigen::VectorXd> const &u0 = std::nullopt,
std::optional<Eigen::VectorXd> const &lambda0 = std::nullopt) {
Eigen::VectorXd u = u0.value_or(Eigen::VectorXd::Zero(rhs.size()));
Eigen::VectorXd u = u0.value_or(0.5*(lowerBound + upperBound));
Eigen::VectorXd lambdas = lambda0.value_or(Eigen::VectorXd::Zero(rhs.size()));
std::vector<bool> is_active_lower(rhs.size(), false);
std::vector<bool> is_active_upper(rhs.size(), false);
std::vector<bool> is_active(rhs.size(), false);
int n_changed = rhs.size() + 20;
int iter = 0;
do {
while(iter < 1000){
//first go through all the conditions to see, which values are active
auto const is_active_old = is_active;
for (size_t i = 0; i < rhs.size(); ++i) {
......@@ -83,9 +83,9 @@ namespace mfds {
is_active_lower[i] = false;
is_active_upper[i] = false;
if (lambdas[i] + alpha * (lowerBound[i] - u[i]) > 0) {
is_active_lower[i] = true;
is_active_lower[i] = !isBoundary[i];
} else if (lambdas[i] + alpha * (upperBound[i] - u[i]) < 0) {
is_active_upper[i] = true;
is_active_upper[i] = !isBoundary[i];
}
is_active[i] = is_active_lower[i] || is_active_upper[i];
}
......@@ -100,8 +100,6 @@ namespace mfds {
}
return n;
}();
std::cout << std::endl;
std::cout << "n_changed: " << n_changed << std::endl;
//rhs adaption has to be done different, as we do not fix values
Eigen::VectorXd cur_rhs = rhs;
{
......@@ -144,11 +142,26 @@ namespace mfds {
lambdas[i] = 0.;
}
}
std::cout << std::endl;
std::cout << "n_changed: " << n_changed << std::endl;
if( (n_changed == 0 && iter > 0) || iter > 100){
if(iter>100){
std::cout << "Max Newton iterations reached" << std::endl;
}
break;
}
//calculate difference to old solution
auto const diff_u = u - u_old;
auto const diff_lambdas = lambdas - lambdas_old;
std::cout << "Diff u: " << diff_u.norm() << ", Diff lambdas: " << diff_lambdas.norm() << std::endl;
} while (n_changed > 0 || iter++ < 1);
iter++;
double const omega = 0.5;
u = u_old + omega*diff_u;
lambdas = lambdas_old + omega*diff_lambdas;
}
return u;
}
......@@ -211,6 +224,10 @@ namespace mfds {
}
}
Eigen::VectorXd const& get_diagonal_preconditioner() const{
return op->get_diagonal_preconditioner();
}
private:
Op_t *op; //just use pointer and leave responsibility to user
std::vector<bool> isBoundary;
......@@ -229,6 +246,8 @@ namespace mfds {
double operator()(size_t i, size_t j) const;
Eigen::VectorXd const& get_diagonal_preconditioner() const;
protected:
void construct_spatial_matrices(mfem::FiniteElementSpace &X);
......@@ -243,6 +262,7 @@ namespace mfds {
Eigen::SparseMatrix<double> mass_t;
std::vector<int> space_time_boundary_dofs;
Eigen::VectorXd generalized_eigenvalues;
Eigen::VectorXd mass_diagonal;
double rho; //regularization parameter
};
......
mfds/conjugate_gradient.hpp
+ 16
− 10
View file @ fbf503e8
......@@ -25,9 +25,12 @@ template<class Op_t, class Vec_t>
class diagonal_preconditioner_t {
public:
diagonal_preconditioner_t(Op_t const &A_, int const n) : A(A_) {
diag = get_diagonal_brute_force(A_, n);
//diag = get_diagonal_brute_force(A_, n);
diag = std::vector<double>(n, 1.);
};
diagonal_preconditioner_t(Op_t const &A_,std::vector<double> const diag_) : A(A_), diag(diag_) {};
Vec_t operator*(Vec_t const &v) const {
Vec_t res = A*v;
for (int i = 0; i < diag.size(); ++i) {
......@@ -56,22 +59,25 @@ Vec_t conjugate_gradient(Op_t const &A_, Vec_t const &b_, Vec_t const &x0, doubl
Vec_t b = A.modify_rhs(b_);
Vec_t x = x0;
Vec_t r = b - A * x;
Vec_t p = r;
Vec_t Ap = A * p;
Vec_t d = r;
Vec_t z = A * d;
double alpha, beta, r_norm, r_norm_new;
r_norm = r.norm();
double const b_norm = b.norm();
std::cout << "b_norm: " << b_norm << std::endl;
for (int i = 0; i < max_iter; ++i) {
alpha = r.dot(r) / p.dot(Ap);
x = x + alpha * p;
r = r - alpha * Ap;
alpha = r.dot(r) / d.dot(z);
x = x + alpha * d;
r = r - alpha * z;
r_norm_new = r.norm();
//std::cout << "Iteration " << i << " Residual: " << r_norm_new << std::endl;
std::cout << "Iteration " << i << " (rel.) Residual: " << r_norm_new << std::endl;
if (r_norm_new < tol) {
break;
}
beta = r.dot(r) / r_norm;
p = r + beta * p;
Ap = A * p;
beta = r_norm_new / r_norm;
beta = beta * beta;
d = r + beta * d;
z = A * d;
r_norm = r_norm_new;
if (i == max_iter - 1) {
std::cerr << "Conjugate gradient did not converge after " << max_iter
......
mfds/tests/active_set_wrapper_test.cpp
+ 12
− 4
View file @ fbf503e8
......@@ -72,16 +72,24 @@ TEST_CASE("Active Set Wrapper", "[mfds][active_set]"){
is_active[7] = true;
mfds::active_set_lin_op_t op(&M, is_active);
SECTION("Check if it does, what is expected on Matrix level"){
Eigen::VectorXd const v = Eigen::VectorXd::Ones(ndof);
Eigen::VectorXd const v = Eigen::VectorXd::Random(ndof)*10;
Eigen::VectorXd const res = op*v;
Eigen::VectorXd const res_mod = get_modified_matrix(M, is_active)*v;
Eigen::VectorXd const diff = res - res_mod;
REQUIRE(diff.norm() == 0.);
REQUIRE(diff.norm() < 1e-8);
if(false){
//print eigenvalues of modified matrix
Eigen::EigenSolver<Eigen::MatrixXd> es(get_modified_matrix(M, is_active));
std::cout << "Eigenvalues of modified matrix: " << std::endl;
std::cout << es.eigenvalues() << std::endl;
}
}SECTION("Check if CG converges"){
Eigen::VectorXd const rhs = Eigen::VectorXd::Random(ndof);
Eigen::VectorXd const u0 = Eigen::VectorXd::Zero(ndof);
Eigen::VectorXd const res = conjugate_gradient(op, rhs, u0, 1e-7, 1000);
Eigen::VectorXd const diff = res - get_modified_matrix(M, is_active).inverse()*rhs;
//Eigen::VectorXd const res = conjugate_gradient(op, rhs, u0, 1e-7, 1000);
Eigen::VectorXd const res = conjugate_gradient(get_modified_matrix(M, is_active), rhs, u0, 1e-7, 1000);
Eigen::VectorXd const res_expected = get_modified_matrix(M, is_active).inverse()*rhs;
Eigen::VectorXd const diff = res - res_expected;
REQUIRE(diff.norm() < 1e-6);
}
}
......
playground/example_matrix_free_operator_2d.cpp
+ 18
− 14
View file @ fbf503e8
......@@ -15,7 +15,8 @@
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);
//return 16*t * t * x * (1 - x) * y * (1 - y);
return std::sin(M_PI*x)*std::sin(M_PI*y)*std::sin(M_PI*t);
}
static double lower_bound_fun(double const t, mfem::Vector const &p) {
......@@ -23,20 +24,21 @@ static double lower_bound_fun(double const t, mfem::Vector const &p) {
}
static double upper_bound_fun(double const t, mfem::Vector const &p) {
return .5* target_function(p,t);
return .5;//std::min(target_function(p, t), .5);
}
int main(int argc, char *argv[]) {
int const nx = 20;
int const ny = 20;
int const nt = 20;
int const refinement = 3;
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-10;
int const cg_max_iter = 100;
double const cg_tol = 1e-12;
int const cg_max_iter = 500;
double const newton_relax_param = 1;
mfem::Mesh spatial_mesh = mfem::Mesh::MakeCartesian2D(nx, nx, mfem::Element::TRIANGLE, true, 1.0, 1.0);
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);
......@@ -53,13 +55,13 @@ int main(int argc, char *argv[]) {
//get space time boundary
auto const space_time_boundary_dofs = heat_operator.get_space_time_boundary_dofs();
std::vector<bool> isBoundary(X.GetTrueVSize() * nt, false);
//std::cout << "Boundary dofs:" << std::endl;
std::cout << "Boundary dofs:" << std::endl;
for (auto const dof: space_time_boundary_dofs) {
//std::cout << dof << std::endl;
std::cout << dof << std::endl;
isBoundary[dof] = true;
}
//calculate inhomogenious rhs for target
//calculate inhomogeneous rhs for target
auto const rhs_inhom = mfds::compute_spatially_inhom_rhs(target_function, T_test, X);
//wrap operator so that boundary conditions can be applied
......@@ -79,10 +81,11 @@ int main(int argc, char *argv[]) {
//solve via active set strategy (keep values consistent)
Eigen::VectorXd u0 = 0.5*(lower_bound + upper_bound);
Eigen::VectorXd u0 = .5*(lower_bound + upper_bound);
Eigen::VectorXd lambda0 = rhs_boundary_adapted - boundary_condition_wrapper*u0;
lambda0 *= 0.;
auto const u_vec = active_set_solver.solve(cg_tol, cg_max_iter, newton_relax_param);
auto const u_vec = active_set_solver.solve(cg_tol, cg_max_iter, newton_relax_param, u0, lambda0);
//add initial condition
auto const state_block = vtint::split_to_block_vector(u_vec, nt);
......@@ -98,12 +101,13 @@ int main(int argc, char *argv[]) {
};
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);
std::cout << "L2 error: " << target_state_l2_error << std::endl;
return 0;
}
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment