Skip to content

Instantly share code, notes, and snippets.

@precise-simulation
Created September 28, 2021 09:14
Show Gist options
  • Save precise-simulation/720cc6ab0517f554684953b42929c15e to your computer and use it in GitHub Desktop.
Save precise-simulation/720cc6ab0517f554684953b42929c15e to your computer and use it in GitHub Desktop.
FEATool solver hook example for implementing the (summation) constraint c1 = c2 over a subdomain
function [ fea ] = featool_solver_hook_example()
%FEATOOL_SOLVER_HOOK_EXAMPLE
%
% Solver hook example to add a constraint that sum(ci_j) = 0 for two
% dependent variables c1 and c2 (and j is the nodes/degrees of freedom).
% Copyright 2021 Precise Simulation, Ltd.
% Geometry and grid.
fea.sdim = { 'x', 'y' };
fea.geom.objects = { gobj_rectangle() };
fea.grid = rectgrid(10);
% Equation settings.
fea = addphys( fea, @convectiondiffusion, { 'c1' } );
fea.phys.cd.eqn.coef = { 'dts_cd', 'd_ts', 'Time scaling coefficient', { '1' };
'd_cd', 'D', 'Diffusion coefficient', { '1' };
'u_cd', 'u', 'Convection velocity in x-direction', { '0' };
'v_cd', 'v', 'Convection velocity in y-direction', { '0' };
'r_cd', 'R', 'Reaction rate', { '1' };
'c10_cd', 'c1_0', 'Initial condition for c1', { '0' } };
fea = addphys( fea, @convectiondiffusion, { 'c2' } );
fea.phys.cd2.eqn.coef = { 'dts_cd2', 'd_ts', 'Time scaling coefficient', { '1' };
'd_cd2', 'D', 'Diffusion coefficient', { '1' };
'u_cd2', 'u', 'Convection velocity in x-direction', { '0' };
'v_cd2', 'v', 'Convection velocity in y-direction', { '0' };
'r_cd2', 'R', 'Reaction rate', { '0' };
'c20_cd2', 'c2_0', 'Initial condition for c2', { '0' } };
% Boundary settings.
fea.phys.cd.bdr.sel = [ 1, 1, 2, 2 ];
fea.phys.cd.bdr.coef = { 'bcr_cd', 'c1 = c1_0', 'Concentration', { 'c1_0' }, { 1, 1, 1, 1 }, [], { '0', '0', '0', '0' };
'bcc_cd', 'n.(-Dgrad c1) = 0', 'Convective flux/outflow', [], { 0, 0, 0, 0 }, [], { 0, 0, 0, 0 };
'bci_cd', 'n.(-Dgrad c1 + uc1) = 0', 'Insulation/symmetry', [], { 0, 0, 0, 0 }, { '(nx*u_cd+ny*v_cd)*c1' }, { '0', '0', '0', '0' };
'bcf_cd', '-n.(-Dgrad c1 + uc1) = N_0', 'Flux condition', { 'N_0' }, { 0, 0, 0, 0 }, { '(nx*u_cd+ny*v_cd)*c1' }, { '0', '0', '0', '0' } };
fea.phys.cd2.bdr.sel = [ 2, 2, 1, 1 ];
fea.phys.cd2.bdr.coef = { 'bcr_cd2', 'c2 = c2_0', 'Concentration', { 'c2_0' }, { 1, 1, 1, 1 }, [], { 0, 0, '1-x', 'y' };
'bcc_cd2', 'n.(-Dgrad c2) = 0', 'Convective flux/outflow', [], { 0, 0, 0, 0 }, [], { 0, 0, 0, 0 };
'bci_cd2', 'n.(-Dgrad c2 + uc2) = 0', 'Insulation/symmetry', [], { 0, 0, 0, 0 }, { '(nx*u_cd2+ny*v_cd2)*c2' }, { '0', '0', '0', '0' };
'bcf_cd2', '-n.(-Dgrad c2 + uc2) = N_0', 'Flux condition', { 'N_0' }, { 0, 0, 0, 0 }, { '(nx*u_cd2+ny*v_cd2)*c2' }, { '0', '0', '0', '0' } };
% Define solver hook to call before/after solving linear system.
fea.sol.settings.hook = @my_solve_hook;
% Solver call.
fea = parsephys( fea );
fea = parseprob( fea );
fea.sol.u = solvestat( fea );
% Postprocessing.
subplot(1,3,1)
postplot( fea, 'surfexpr', 'c1', 'title', 'Concentration, c1' );
subplot(1,3,2)
postplot( fea, 'surfexpr', 'c2', 'title', 'Concentration: c2' );
subplot(1,3,3)
postplot( fea, 'surfexpr', 'c1+c2', 'title', 'c1+c2' );
% Check results.
ind_c1 = find(strcmp(fea.dvar, 'c1'));
offset = sum(fea.eqn.ndof(1:ind_c1-1));
ind_dof_c1 = unique(fea.eqn.dofm{ind_c1}(:)) + offset;
ind_c2 = find(strcmp(fea.dvar, 'c2'));
offset = sum(fea.eqn.ndof(1:ind_c2-1));
ind_dof_c2 = unique(fea.eqn.dofm{ind_c2}(:)) + offset;
c1_plus_c2 = sum(fea.sol.u(ind_dof_c1)) + ...
sum(fea.sol.u(ind_dof_c2))
assert( abs(c1_plus_c2) < 1e-12 )
%------------------------------------------------------------------------------%
function [ A, f, prob ] = my_solve_hook( A, f, prob, solve_step )
%MY_SOLVE_HOOK Custom solver hook calling function.
%
% [ A, F, PROB ] = SOLVE_HOOK( A, F, PROB, SOLVE_STEP )
% Special solver hook function. Allows for special functions to be
% called before and after the linear solver to modify the global
% matrix A and right hand side/load vector F with the information in
% the finite element problem struct PROB. SOLVE_STEP indicates if
% the function call is before the linear solver call (=1), or to
% clean-up after (=2).
%
% A solve hook may be prescribed by entering a function name string
% or handle in the fea.sol.settings.hook field. The corresponding
% function will in each solve iteration then be called as
%
% [A,f,fea] = fcn( A, f, fea, solve_step );
if( solve_step == 1 )
disp('Calling solver hook BEFORE linear solver call.')
ind_c1 = find(strcmp(prob.dvar, 'c1'));
offset = sum(prob.eqn.ndof(1:ind_c1-1));
ind_dof_c1 = unique(prob.eqn.dofm{ind_c1}(:)) + offset;
n = length(ind_dof_c1);
ind_c2 = find(strcmp(prob.dvar, 'c2'));
offset = sum(prob.eqn.ndof(1:ind_c2-1));
ind_dof_c2 = unique(prob.eqn.dofm{ind_c2}(:)) + offset;
if( length(ind_dof_c1) ~= length(ind_dof_c2) )
error('c1 and c2 must have the same discretization (fea.sfun)')
end
if( isstruct(A) ) % Sparsify matrix if necessary.
A = sparse( A.rows, A.cols, A.vals, A.size(1), A.size(2) );
end
% Build constraint contributions.
B = sparse(1, size(A,2));
B(sub2ind(size(B), ones(1,n), ind_dof_c1(:)')) = 1;
B(sub2ind(size(B), ones(1,n), ind_dof_c2(:)')) = 1;
% Add row to add constraint sum_i(c1_i + c2_i) = 0.
A = [A, B';
B, 0 ];
f = [f; 0];
end
if( solve_step == 2 )
disp('Calling solver hook AFTER linear solver call.')
% Remove system modifications.
n0 = sum(prob.eqn.ndof);
A = A(1:n0,1:n0);
f = f(1:n0);
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment