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

5 messages
Open this post in threaded view
|

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

 This post was updated on . 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.
Open this post in threaded view
|

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

 Administrator 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.
Open this post in threaded view
|

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

 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