From 5dc132a6f33617d7858b8c1b5d5ab5b1d2e2215e Mon Sep 17 00:00:00 2001 From: Rafael Vazquez Date: Mon, 22 Jul 2024 13:39:24 +0200 Subject: [PATCH] Moved functions to separate files --- .../op_nitsche_consistency_cahn_hilliard.m | 2 +- .../multipatch/mp_solve_cahn_hilliard_C1.m | 162 +++++++++--------- .../solve/cahn_hilliard/Res_K_cahn_hilliard.m | 43 +++++ .../generalized_alpha_step_cahn_hilliard.m | 95 ++++++++++ geopdes/inst/solve/solve_cahn_hilliard.m | 120 ++++++------- .../op_nitsche_consistency_cahn_hilliard.m | 2 +- 6 files changed, 281 insertions(+), 143 deletions(-) create mode 100644 geopdes/inst/solve/cahn_hilliard/Res_K_cahn_hilliard.m create mode 100644 geopdes/inst/solve/cahn_hilliard/generalized_alpha_step_cahn_hilliard.m diff --git a/geopdes/inst/multipatch/@sp_multipatch_C1/op_nitsche_consistency_cahn_hilliard.m b/geopdes/inst/multipatch/@sp_multipatch_C1/op_nitsche_consistency_cahn_hilliard.m index f67be7b9..ec4e918c 100644 --- a/geopdes/inst/multipatch/@sp_multipatch_C1/op_nitsche_consistency_cahn_hilliard.m +++ b/geopdes/inst/multipatch/@sp_multipatch_C1/op_nitsche_consistency_cahn_hilliard.m @@ -29,7 +29,7 @@ % You should have received a copy of the GNU General Public License % along with this program. If not, see . -function [A] = op_nitsche_consistency_cahn_hilliard (space, msh, nmnn_sides, coeff) +function A = op_nitsche_consistency_cahn_hilliard (space, msh, nmnn_sides, coeff) if (nargin == 3 || isempty(coeff)) coeff = @(varargin) ones (size(varargin{1})); diff --git a/geopdes/inst/multipatch/mp_solve_cahn_hilliard_C1.m b/geopdes/inst/multipatch/mp_solve_cahn_hilliard_C1.m index 3dc75453..b9a32932 100644 --- a/geopdes/inst/multipatch/mp_solve_cahn_hilliard_C1.m +++ b/geopdes/inst/multipatch/mp_solve_cahn_hilliard_C1.m @@ -168,7 +168,7 @@ disp('----------------------------------------------------------------') disp(strcat('time step t=',num2str(time))) - [u_n1, udot_n1] = generalized_alpha_step(u_n, udot_n, dt, a_m, a_f, gamma, mu, dmu, ... + [u_n1, udot_n1] = generalized_alpha_step_cahn_hilliard (u_n, udot_n, dt, a_m, a_f, gamma, mu, dmu, ... mass_mat, lapl_mat, bnd_mat, Pen, pen_rhs, space, msh); % check flux through the boundary @@ -215,88 +215,88 @@ %-------------------------------------------------------------------------- +% %-------------------------------------------------------------------------- +% % One step of generalized alpha-method +% %-------------------------------------------------------------------------- +% +% function [u_n1, udot_n1] = generalized_alpha_step(u_n, udot_n, dt, a_m, a_f, gamma, mu, dmu, ... +% mass_mat, lapl_mat, bnd_mat, Pen, pen_rhs, space, msh) +% +% % Convergence criteria +% n_max_iter = 20; +% tol_rel_res = 1e-10; +% tol_abs_res = 1e-10; +% +% % Predictor step +% u_n1 = u_n; +% udot_n1 = (gamma-1)/gamma * udot_n; +% +% % Newton loop +% for iter = 0:n_max_iter +% +% % Field at alpha level +% udot_a = udot_n + a_m *(udot_n1-udot_n); +% u_a = u_n + a_f *(u_n1-u_n); +% +% % Compute the residual (internal) +% [Res_gl, stiff_mat] = Res_K_cahn_hilliard (space, msh, ... +% mass_mat, lapl_mat, bnd_mat, Pen, pen_rhs, ... +% u_a, udot_a, mu, dmu); +% +% % Convergence check +% if (iter == 0) +% norm_res_0 = norm(Res_gl); +% end +% norm_res = norm(Res_gl); +% +% if (norm_res/norm_res_0 < tol_rel_res) +% disp(strcat('iteration n°=',num2str(iter))) +% disp(strcat('norm (abs) residual=',num2str(norm_res))) +% break +% end +% if (norm_res. + +function [Res_gl, stiff_mat] = Res_K_cahn_hilliard (space, msh, mass_mat, lapl_mat, ... + bnd_mat, Pen, pen_rhs, u_a, udot_a, mu, dmu) + + % Double well (matrices) + if (isa (space, 'sp_scalar')) + [term2, term2K] = op_gradfu_gradv_tp (space, msh, u_a, mu, dmu); + elseif (isa (space, 'sp_multipatch_C1')) + [term2, term2K] = op_gradfu_gradv_mp (space, msh, u_a, mu, dmu); + else + error ('Not implemented for this kind of space') + end + + % Residual + Res_gl = mass_mat*udot_a + term2*u_a + lapl_mat*u_a; + + % Tangent stiffness matrix (mass is not considered here) + stiff_mat = term2 + term2K + lapl_mat; + + % In case of Neumann BC, add boundary terms + if (~isempty(bnd_mat)) + Res_gl = Res_gl - (bnd_mat + bnd_mat.') * u_a + Pen*u_a - pen_rhs; + stiff_mat = stiff_mat - (bnd_mat + bnd_mat.') + Pen; + end +end diff --git a/geopdes/inst/solve/cahn_hilliard/generalized_alpha_step_cahn_hilliard.m b/geopdes/inst/solve/cahn_hilliard/generalized_alpha_step_cahn_hilliard.m new file mode 100644 index 00000000..08462666 --- /dev/null +++ b/geopdes/inst/solve/cahn_hilliard/generalized_alpha_step_cahn_hilliard.m @@ -0,0 +1,95 @@ +% GENERALIZED_ALPHA_STEP_CAHN_HILLIARD: perform one step of the generalized alpha method +% for the solution of the Cahn-Hilliard equation, solving the nonlinear equation with Newton's method. +% It is called from solve_cahn_hilliard and mp_solve_cahn_hilliard_C1. +% +% INPUT: +% +% u_n: field at the previous time step. +% u_dotn: time derivative at the previous time step. +% dt: time step size +% a_m: parameter for the generalized alpha method +% a_f: parameter for the generalized alpha method +% gamma: parameter for the generalized alpha method +% mu, dmu: function handle for value and derivative of the mu coefficient +% mass_mat: mass matrix, to avoid recomputing +% lapl_mat: bilaplacian matrix (Delta u, Delta v), to avoid recomputing +% bnd_mat: matrix with Nitsche's consistency term, to avoid recomputing +% Pen: penalty matrix from Nitsche's method +% pen_rhs: penalty rhs from Nitsche's method +% space: space object (see sp_scalar or sp_multipatch_C1) +% msh: mesh object (see msh_cartesian or msh_multipatch) +% +% OUTPUT: +% +% u_n1: field at the new step. +% u_dotn1: time derivative at the new step. +% +% Copyright (C) 2023, 2024 Michele Torre, Rafael Vazquez +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. + +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + +function [u_n1, udot_n1] = generalized_alpha_step_cahn_hilliard(u_n, udot_n, dt, a_m, a_f, gamma, mu, dmu, ... + mass_mat, lapl_mat, bnd_mat, Pen, pen_rhs, space, msh) + +% Convergence criteria + n_max_iter = 20; + tol_rel_res = 1e-10; + tol_abs_res = 1e-10; + +% Predictor step + u_n1 = u_n; + udot_n1 = (gamma-1)/gamma * udot_n; + +% Newton loop + for iter = 0:n_max_iter + + % Field at alpha level + udot_a = udot_n + a_m *(udot_n1-udot_n); + u_a = u_n + a_f *(u_n1-u_n); + + % Compute the residual (internal) + [Res_gl, stiff_mat] = Res_K_cahn_hilliard (space, msh, ... + mass_mat, lapl_mat, bnd_mat, Pen, pen_rhs, ... + u_a, udot_a, mu, dmu); + + % Convergence check + if (iter == 0) + norm_res_0 = norm(Res_gl); + end + norm_res = norm(Res_gl); + + if (norm_res/norm_res_0 < tol_rel_res) + disp(strcat('iteration n°=',num2str(iter))) + disp(strcat('norm (abs) residual=',num2str(norm_res))) + break + end + if (norm_res