Why Periodic BCs does not converge when using the Newton solver?

Hello，
   I encountered difficulties in using Periodic BCs. In two dependent variable case, Periodic BCs  can only converge in very simple equations (one time Newton to converge). If the equations are nonlinear, it does not converge. For example, in a rectangle(0, 1, 0, 1) region, two simple Poisson- equations. It can converge with all DirBCs. If changing to Periodic BCs (Left =Right, using ex_swirl_flow3.m), it does not converge. The code is


fea.phys.ce.eqn.seqn = 'u'' + (ux_x + uy_y) = sin(x)+q';
fea.phys.ce2.eqn.seqn = 'q'' + (qx_x + qy_y)=sin(x)*u';
fea.phys.ce.bdr.coef = { 'bcnd_ce', '', '', {''}, ...
                         { 1, 0,      1, 0 }, [], ...
                         { 0    , @periodic_bc, 0,    0     } };

fea.phys.ce2.bdr.coef = { 'bcnd_ce2', '', '', {''}, ...
                         { 1, 0,      1, 0 }, [], ...
                         { 0    , @periodic_bc, 0,    0    } };

The Newton solver:

fea.sol.u = solvestat( fea, 'nsolve',2);


if( i_bdr~=2 || solve_step~=1 )  
  return                        
end
j_bdr = 4;

If 'nsolve' to be 1 or -1, it can converge. However, if using the Newton solver, it does not converge. I can not find the reason. The pure Newton solver is important for other related problems.

It is not really clear what example you are referring to. swirl_flow3 does include periodic boundary conditions on top and bottom but is a time dependent problem and does not use the Newton method. If possible, please include a *full* example to diagnose issues.

close all
clc
clear
fea.sdim = {'x' 'y'};
gobj = gobj_rectangle( 0, 1, 0, 1, 'R1' );
fea = geom_add_gobj( fea, gobj );
fea.grid = gridgen( fea, 'gridgen', 'default','quad',0, 'hmax', 0.1, 'grading', 0.3 );

%% Equation settings
fea = addphys( fea, @customeqn, { 'u' } );
fea.phys.ce.dvar = { 'u' };
fea.phys.ce.eqn.seqn = 'u'' + (ux_x + uy_y) = sin(x)+q';

fea = addphys( fea, @customeqn, { 'q' } );
fea.phys.ce2.dvar = { 'q' };
fea.phys.ce2.eqn.seqn = 'q'' + (qx_x + qy_y)=sin(x)*u';

fea.phys.ce.bdr.coef = { 'bcnd_ce', '', '', {''}, ...
                         { 1, 0,      1, 0 }, [], ...
                         { 0    , @periodic_bc2, 0,    0     } };


fea.phys.ce2.bdr.coef = { 'bcnd_ce2', '', '', {''}, ...
                         { 1, 0,      1, 0 }, [], ...
                         { 0    , @periodic_bc2, 0,    0    } };


% fea.phys.ce.bdr.coef = { 'bcnd_ce', '', '', {''}, ...
%                          { 1, 1,      1, 1 }, [], ...
%                          { 0    , 0, 0,    0     } };
%
%
% fea.phys.ce2.bdr.coef = { 'bcnd_ce2', '', '', {''}, ...
%                          { 1, 1,      1, 1 }, [], ...
%                          { 0    , 0, 0,    0    } };


fea = parsephys( fea );
fea = parseprob( fea );
fea.sol.u = solvestat( fea, 'nsolve',1);
% fea.sol.u = solvestat( fea, 'nsolve',2);


figure;
postplot( fea, 'surfexpr', 'u', 'title', 'surface: Dependent variable, u', 'solnum', 1 );

figure
clf
postplot( fea, 'surfexpr', 'q', 'title', 'surface: Dependent variable, q','solnum', 1 );













function [ A, f, prob ] = periodic_bc2( A, f, prob, i_dvar, i_bdr, solve_step )
if( i_bdr~=2 || solve_step~=1 )   % Only process boundaries 1 ( and 3 ) in
  return                          % the step directly before linear solver.
end


j_bdr = 4;

if( isstruct(A) )
  A = sparse( A.rows, A.cols, A.vals, A.size(1), A.size(2) );
end

bdrm = prob.bdr.bdrm{i_dvar};
ndof = prob.eqn.ndof;

dof_offset = sum( prob.eqn.ndof(1:(i_dvar-1)) );


ix_i     = find( bdrm(3,:)==i_bdr );   % Index to dofs on boundary i_bdr.
[~,itmp] = unique(bdrm(4,ix_i));       % Remove duplicate/shared points (optional).
ix_i     = sort( ix_i(itmp) );         % Sort index.
for i=1:numel(ix_i)
  ix = ix_i(i);
  x_i(i) = evalexpr0( 'y', bdrm(6:end,ix), [], bdrm(1,ix), bdrm(2,ix), prob );
end


ix_j     = find( bdrm(3,:)==j_bdr );   % Index to dofs on boundary j_bdr.
[~,itmp] = unique(bdrm(4,ix_j));       % Remove duplicate/shared points (optional).
ix_j     = sort( ix_j(itmp) );         % Sort index.
ix_ij    = zeros( size(ix_i) );        % Index linking i and j dofs.
for j=1:numel(ix_j)
  ix = ix_j(j);
  x_j = evalexpr0( 'y', bdrm(6:end,ix), [], bdrm(1,ix), bdrm(2,ix), prob );

  [~,ix_ij_j] = min( abs( x_i - x_j ) );
  ix_ij(ix_ij_j) = ix;
end
idofs = bdrm(4,ix_i)  + dof_offset;


jdofs = bdrm(4,ix_ij) + dof_offset;


[~,imin] = min(x_i);
[~,imax] = max(x_i);
ind_rem = [imin imax];
idofs(ind_rem) = [];
jdofs(ind_rem) = [];

A(jdofs,:) = A(idofs,:) + A(jdofs,:);
f(jdofs)   = f(idofs)   + f(jdofs);

A(idofs,:) = 0;
A(sub2ind(size(A),idofs,idofs)) =  1;
A(sub2ind(size(A),idofs,jdofs)) = -1;

f(idofs) = 0;



end

Thank you for reporting. We have identified an issue with the Newton solver and solver hooks, and will include a correction in the next update later this year.

Thank you. Looking forward to the next release.