residual.m

来自「MDPSAS工具箱是马里兰大学开发的」· M 代码 · 共 115 行

function [A, Jacobian, Rescol] = residual(varargin)

% residual.m Nonlinear equation model object residual update.
%
% This method is for storing implicit equations and updating 'resid'
% data field for nonlinear eqn models, for example, f(x,y) = 0
% 
% If classes derived from solver do not include an evalvar method,
% this method should be overloaded.

A = varargin{1};

if length(varargin) == 7
   Akey       = varargin{2};
   Aval       = varargin{3};
   V          = varargin{4};
   template   = varargin{5};
   Epsil      = varargin{6};
   VPlusEpsil = varargin{7};
   end

unpack(A,'var')
unpack(A,'param')

if length(varargin) == 7 % update residual and compute jacobian
   ntot = length(V);
else                     % only update residual
   ntot = 0;
   end

i = 1;  % key index
j = 1;  % value index

for n = 0:ntot

   if n > 0
      if isa(eval(Akey{i}),'double')
         eval([Akey{i},'(j) = VPlusEpsil(n);']);
      elseif isa(eval(Akey{i}),'scalarfield')
         eval([Akey{i},' = setval(',Akey{i},'VPlusEpsil(n),j);']);
      else
         error('Unsupported variable class')
         end
      end

% .......... cut and replace everything below .......... %

vx    = a*Bvx;
vy    = b*Bvy + vyBC;
press = d*Bp;

% Interior residuals

Rescell{1} = wip( DDx*vx + DDy*vx - Dx*press, Bvx );
Rescell{2} = wip( DDx*vy + DDy*vy - Dy*press, Bvy );
Rp  = wip( Dx*vx + Dy*vy, Bp );
Rp(1,1) = d(1,1)-0; % set mean pressure = 0;
Rescell{3} = Rp;

% No boundary condition residuals

% .......... cut and replace everything above .......... %

   if n == 0
      A = set(A,'resid',Rescell);
      Rescol = residcol(Rescell);
      end

  if n > 0
      if isa(eval(Akey{i}),'double')
         eval([Akey{i},'(j) = V(n);']);
      elseif isa(eval(Akey{i}),'scalarfield')
         eval([Akey{i},' = setval(',Akey{i},'V(n),j);']);
      else
         error('Unsupported variable class')
         end

      j = j+1;
      if j > prod(template{i})  % move to the next key
         j = 1;
         i = i+1;
         end

      RPlusEpsil = residcol(Rescell);
      Jacobian(:,n) = ( RPlusEpsil - Rescol )/Epsil(n);
      end

   end

% -------------------- residcol.m -------------------- %

function column = residcol(valuecell)

% Convert a cell vector (row or column format) of doubles and/or
% scalarfields to a column.

noel = length(valuecell);
column = [];

for i = 1:noel

   if isa(valuecell{i},'double')
      temp = size(valuecell{i});
   elseif isa(valuecell{i},'scalarfield')
      temp = valsize(valuecell{i});
    else
      error('Must be a double or scalarfield object')
      end

   column = [column; reshape(valuecell{i},prod(temp),1)];
   end