ANSYS FLUENT 12.0 Theory Guide - 18.6.4 Full-Approximation Storage (FAS) Multigrid (enea.it)
多重网格时间规约算法（MGRIT）允许对时间演化微分方程的求解进行时间域并行化。该过程基于多重网格规约技术：更多信息可以在原始MGRIT论文中找到： [1] "Parallel time integration with multigrid", R. Falgout, S. Friedhoff, Tz. V. Kolev, S. P. MacLachlan, and J. B. Schroder (https://computation.llnl.gov/projects/parallel-time-integration-multigrid/mgritPaper-2013-3.pdf)
% In short, the technique is rooted in multigrid reduction: the temporal
% nodes are divided in "Fine" and "Coarse":
% C F F F C F F F C F F F C
% o---x---x---x---o---x---x---x---o---x---x---x---o---> t
% t0 t1 t2 t3 t4 ... tf时间节点分为“细”和“粗”：两个连续 C 节点（例如，[t0, t4]）之间的跨度定义了一个时间块。给定C节点的解，可以通过时间步进（F 松弛）恢复内部 F 节点的解，这个过程对于每个时间块是独立的。这就是时间并行化的潜力所在，因为不同的处理器可以并行地对不同的时间块进行积分。
C 节点可以再次分为细和粗，做到真正的多重网格过程。在最粗层次上，解通过顺序应用时间步进程序（可能与较细层次使用的不同，在更大的时间步上定义）恢复。层次之间的对话通过逐点注入（限制）和注入+F松弛（插值）完成。非线性方程还实现了 FAS 校正方案。在所有层次上，可以应用F或FCF松弛。整个过程反复进行。
此代码仅限于均匀时间网格，但仍然相当优化且灵活。特别是：
- 已收敛的时间块不会被重新计算。（每次应用 F 松弛时，一个更多的时间块达到收敛。算法跟踪这一点并避免冗余计算）。请注意，规定的最大迭代次数不应超过每个 cycle 应用的F松弛次数（算法在这方面不检查收敛性）。
- 可以规定任何可选择的多重网格 cycle。
- 可以在不同层次上规定不同的时间步进器。
%% Initialisation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
verbose = false;
% Initialise multigrid cycle info:
maxLvl = length(NT);
% - here we generate a list of numbers identifying at which order the levels need to be visited,
% and which action needs to be performed. The code is:
% - i->i : smooth -> ex: 1-s-1
% - i->i+1: restrict -> ex: 1-r-2
% - i->i-1: correct -> ex: 2-c-1
% - i=maxLvl: direct solve
% that is, I smooth (s) if I stay at the same level, restrict (r) if I move to coarser one,
% correct (c) if I move to finer one
switch cycleType.name
case 'V-cycle'
% Descend all the way to the bottom, smoothing at each lvl, climb all the way back up, smoothing at each lvl
% 1-s-1-r-2-s-2-r-3-s-3-r-4-c-3-s-3-c-2-s-1-s-1
% NB: notice that if more smoothing iterations are needed, this can be easily adjusted
levels = [repelem( 1:maxLvl-1, cycleType.it(1)+1 ), maxLvl, repelem( maxLvl-1:-1:1, cycleType.it(2)+1 )];
case 'F-cycle'
% Drop to bottom and slowly climb your way to the top, performing a complete V-cycle at each level:
% 1-r-2-r-3-r-4-s-4-c-3-s-3-r-4-s-4-c-3-s-3-c-2-s-2...
levels = 1:maxLvl;
for i=(maxLvl-1):-1:1
% the repelem defines the number of smoothing iterations at each level; the repmat defines the number of complete V-cycles
temp = [ repelem( i:maxLvl-1, cycleType.it(1)+1 ), maxLvl, repelem( maxLvl-1:-1:i, cycleType.it(2)+1 ) ]; % this is a V-cycle
temp = repmat( temp(1:end-1), [ 1, cycleType.it(3) ] ); % this is a V-cycle repeated cycleType.it(3) times
levels = [levels, temp, i ];
end
case 'Custom'
levels = cycleType.lvls;
if max(levels)~=maxLvl
error('MGRIT: there is a mismatch between the number of levels in the cycle prescribed and the discretisation info provided per level');
end
otherwise
error('MGRIT: Type of MG cycle specified not recognised/implemented');
end
% plot(-levels,'x-'); % to make sure we're doing everything right
% Initialise discretisation-related stuff
t = linspace(T(1),T(2),NT(1)+1); % actual discretisation
M = NT ./ NT([2:end,end]); % ratio fine/coarse node (also, M-1 is the number of fine nodes in a chunk per level)
dT = (T(2)-T(1)) ./ NT(1:end); % time step for solver
NX = size(u0,1); % number of variables in the system
% Assign easier names for solvers
timeStep = cell(maxLvl,1);
for lvl = 1:maxLvl % the min is just to deal with the case length(Solver)==1
% timeStep{lvl} = @( u0, p, i ) Solver{ min(lvl, length(Solver)) } ( u0, Problem, dT(lvl), 1, T(1) + (p-1)*DT(lvl) + (i-1)*dT(lvl) );
% timeStep{lvl} = @( u0, i ) Solver{ min(lvl, length(Solver)) } ( u0, Problem, dT(lvl), 1, T(1) + (i-1)*dT(lvl) );
timeStep{lvl} = @( u0, i ) Solver{ min(lvl, length(Solver)) } ( u0, Problem{min(lvl, length(Problem))}, dT(lvl), 1, T(1) + (i-1)*dT(lvl) );
end
% Compute serial solution if needed
if( sSol )
err = zeros(maxIT+1,NT(2)+1); % initialise error as well
uS = Solver{1}( u0, Problem{1}, dT(1), NT(1), T(1) );
end
% Initialise useful info at all levels:
rhs = cell(maxLvl,1); % right-hand side of system
u = cell(maxLvl,1); % complete solution
cIdx = cell(maxLvl,1); % indices corresponding to (unconverged) coarse nodes
fIdx = cell(maxLvl,1); % indices corresponding to (unconverged) fine nodes
for lvl = 1:maxLvl
rhs{lvl} = zeros( NX, NT(lvl)+1 );
u{lvl} = zeros( NX, NT(lvl)+1 );
cIdx{lvl} = 1 : M(lvl) : (NT(lvl)+1);
fIdx{lvl} = 1 : (NT(lvl)+1);
fIdx{lvl}(cIdx{lvl}) = [];
cIdx{lvl}(1) = []; % initial condition is already known at all levels
end
u{1}(:,1) = u0;
rhs{1}(:,1) = u0;
% Initialise solution using a coarse sweep at second level
% u{1}(:,[1,cIdx{1}]) = Solver{2}( u0, Problem{min(2, length(Problem))}, dT(2), NT(2), T(1) );
% notice that rhs at finest level is always considered to be [u0,0,0,...]'
% notice that fIdx at coarsest level is useless (no distinction btw coarse and fine nodes)
% Initialise output
outIdx = 1;
uOut = zeros( NX, NT(1)+1, length(printIt) );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Main MG cycle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
disp('Beginning Multigrid loop')
end
for MGcycle=1:maxIT
% store error
if( sSol )
err(MGcycle,:) = max( abs( u{1}(:,1:M(1):end)-uS(:,1:M(1):end) ), [], 1 );
end
for step = 1:length(levels)-1
lvl = levels(step);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Base case
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if lvl == maxLvl
if verbose
disp(strcat('- "Exact" solve, level: ', int2str(lvl), ' (bottom)'))
end
% Just perform a coarse sweep to solve AC*uC = bC, with AC = spdiag([-G,I],[-1,0])
% - notice each equation is given by -G(u_{p-1}) + u_{p} = rhs_{p}
% - so, to recover u_{p} I need to first propagate G(u_{p-1}) and then add the rhs to it
for p = cIdx{lvl}
uG = timeStep{lvl}( u{lvl}( :, p-1 ), p-1 ); % propagation of u
u{lvl}( :, p ) = uG( :, end ) + rhs{lvl}( :, p ); % correction
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MG recursion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%----------------------------------------------------------------------
% Smoothen
%----------------------------------------------------------------------
if( levels(step) == levels(step+1) )
if verbose
msg = '- F';
if cycleType.FCF
msg = strcat(msg,'-CF');
end
disp(strcat(msg, ' smoothing, level ', int2str(lvl)))
end
% F - smoothing: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - Integrate over all fine nodes starting from values at coarse nodes
% - This process can be carried in parallel over the time chunks
% - Notice the integration stops right before the following coarse node
for i = fIdx{lvl}
% advance by one time step
uF = timeStep{lvl}( u{lvl}( :, i-1 ), i-1 );
% store solution
u{lvl}( :, i ) = uF( :, end ) + rhs{lvl}( :, i );
end
% remove from the list of fine nodes to update all those coming before the first non-converged coarse
% that is, remove the first chunk
fIdx{lvl} = fIdx{lvl}( fIdx{lvl} > min(cIdx{lvl}) );
% If requested, add also CF - smoothing
if cycleType.FCF
% C - smoothing: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - Integrate over the last step in the chunk to update value at following coarse node
% - This process can be carried in parallel over the time chunks
for p = cIdx{lvl}
% advance by last time step
uF = timeStep{lvl}( u{lvl}( :, p-1 ), p-1 );
% store solution
u{lvl}( :, p ) = uF( :, end ) + rhs{lvl}( :, p );
end
% At this point yet another chunk has reached convergence at this level: there is no need to modify this further
% - Remove the node from the list of coarse nodes to update
cIdx{lvl} = cIdx{lvl}( 2 : end );
% - If global convergence is reached, then there's no need for additional refinement below this level
if isempty(cIdx{lvl})
% check at which step I'll be interpolating to a higher level
nextStep = step + find( levels( (step+1) : end ) < lvl, 1 );
if isempty( nextStep )
% if none is found, then I'm at top level and I've reached convergence: stop the whole MGRIT iterations
MGcycle = maxIT;
else
% else, skip to that interpolation step. Notice I need to backtrack twice because:
% - at the end of this loop, the index is increased
% - interpolation is done *before* moving to a higher level
step = nextStep - 2;
end
break;
end
% F - smoothing: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% - Integrate over the following chunks (so, ignoring the first)
% - This process can be carried in parallel over the time chunks
for i = fIdx{lvl}
% advance by one time step
uF = timeStep{lvl}( u{lvl}( :, i-1 ), i-1 );
% store solution
u{lvl}( :, i ) = uF( :, end ) + rhs{lvl}( :, i );
end
% remove M-1 more nodes, ie (almost) one chunk, from the list of fine nodes to update
fIdx{lvl} = fIdx{lvl}( M(lvl) : end );
end
%----------------------------------------------------------------------
% Restrict
%----------------------------------------------------------------------
elseif( levels(step) < levels(step+1) )
if verbose
disp(strcat('- restricting, level: ', int2str(lvl), '->', int2str(lvl+1)))
end
% Compute residual and FAS correction and inject them to coarser level.
% - In order to compute the residual, the final step of fine integration needs to be performed
% - In order to compute the FAS correction, a coarse integration needs to be performed
% - If this comes after a smoothing at the finest level, then yet another coarse node has reached convergence
% Restricted at coarse level, equations satisfy: AFuF = bF -> - F(uC_{p-1}) + uC_{p} = R(bF)_p
% where uC = R(uF) is the restriction of the fine solution to the coarse level
% The residual is then: res = bF-AFuF = R(bF)_p + F(uC_{p-1}) - uC_{p}
% FAS correction requests for an additional term to be added: AC(uC) -> -G(uC_{p-1}) + uC_{p}
% So at the end of the day, the rhs for the coarse level is given by:
% bC = R(bF) + F(uC) -uC + uC - G(uC) = R(bF) + F(uC) - G(uC)
% - initialise temp variables to store F(uC) and G(uC)
uF0 = u{lvl}( :, cIdx{lvl} );
uG0 = u{lvl}( :, cIdx{lvl} );
% - fill them
for pp = 1:length(cIdx{lvl}) % for each (non-converged) chunk
p = cIdx{lvl}(pp);
prevP = p-M(lvl);
% advance by last time step
uF = timeStep{lvl }( u{lvl}( :, p-1 ), p-1 );
% coarse integrate over this chunk
uG = timeStep{lvl+1}( u{lvl}( :, prevP ), prevP );
% store solution
uF0( :, pp ) = uF( :, end );
uG0( :, pp ) = uG( :, end );
end
% If this restriction was preceded by a relaxation, then at this point yet another chunk has reached
% convergence at this level: there is no need to modify this further
if ( step > 1 ) && ( levels(step-1) == levels(step) )
% - Substitute exact solution in this coarse node
u{lvl}( :, cIdx{lvl}(1) ) = uF0(:,1) + rhs{lvl}( :, cIdx{lvl}(1) );
% - Discard first results from integration (now useless)
uF0 = uF0( :, 2:end );
uG0 = uG0( :, 2:end );
% - Remove the node from the list of coarse nodes to update
cIdx{lvl} = cIdx{lvl}( 2 : end );
end
% - If global convergence is reached, then there's no need for additional refinement below this level
if isempty(cIdx{lvl})
% check at which step I'll be interpolating to a higher level
nextStep = step + find( levels( (step+1) : end ) < lvl, 1 );
if isempty( nextStep )
% if none is found, then I'm at top level and I've reached convergence: stop the whole MGRIT iterations
MGcycle = maxIT;
else
% else, skip to that interpolation step. Notice I need to backtrack twice because:
% - at the end of this loop, the index is increased
% - interpolation is done *before* moving to a higher level
step = nextStep - 2;
end
break;
end
% - update converged nodes at sublevels, too
for subLvl = lvl+1:maxLvl
fIdx{subLvl} = ( cIdx{subLvl-1} - 1 ) / M(subLvl-1) + 1; % collect all remaining nodes idx at sublvl
temp = length(fIdx{subLvl}):-M(subLvl):1; % of these, the coarse ones are intervalled as such
cIdx{subLvl} = fIdx{subLvl}( temp(end:-1:1) ); % extract the coarse ones (in ascending order)
fIdx{subLvl}( temp ) = []; % exclude coarse nodes from fine ones
end
% - finally, compute and inject residual + FAS correction
rhs{lvl+1}( :, ((cIdx{lvl}(1)-1)/M(lvl)+1) : end ) = rhs{lvl}( :, cIdx{lvl} ) + uF0 - uG0;
% - also initialise correction at coarser level
% - this is *not* equivalent to taking an additional C-relaxation,
% since uG0 had been computed with the *previous* level at the
% coarse nodes; regardless, it might improve the initial guess
% on the correction at the coarser level:
u{lvl}(:,cIdx{lvl}) = uF0 + rhs{lvl}(:,cIdx{lvl});
% - otherwise comment above to just take the current values of u (including the -eventually- freshly converged one)
u{lvl+1}( :, max((cIdx{lvl}(1)-1)/M(lvl),1) : end ) = u{lvl}( :, [cIdx{lvl}-M(lvl),end] );
%----------------------------------------------------------------------
% Correct
%----------------------------------------------------------------------
elseif( levels(step) > levels(step+1) )
if verbose
disp(strcat('- interpolating, level: ', int2str(lvl), '->', int2str(lvl-1)))
end
% Correction is made via ideal interpolation:
% - Collect the values of the solution at this level
% - Integrate at the finer level starting from these values (this is equivalent to F-relaxation at higher level)
% - This makes the complete correction: *substitute* it to the solution at the finer level
% - Notice the rhs *should* appear
% Correct coarse points first: since injection is used, things simplify a bit, and I can just inject back
u{lvl-1}(:,cIdx{lvl-1}) = u{lvl}( :, (end-length(cIdx{lvl-1})+1):end);
% If there is a smoothing afterwards, using ideal interpolation is useless (its result gets overwritten),
% so only do it if that's not the case
if ( step == length(levels)-1 ) || ( levels(step+1) ~= levels(step+2) )
for i = fIdx{lvl-1} % for each chunk at higher level (ie, each node at this level)
% advance by one time step
uF = timeStep{lvl-1}( u{lvl-1}(:,i-1), i-1 );
% store solution
u{lvl-1}( :, i ) = uF( :, end ) + rhs{lvl-1}( :, i );
end
end
end
% end of MG cycle
end
if verbose
disp(strcat('End of MG cycle #', int2str(MGcycle)))
end
% print temporary results
if(MGcycle == printIt(outIdx) )
uOut(:,:,outIdx) = u{1};
outIdx = min(outIdx+1, length(printIt));
end
% if(~mod(MGcycle,10))
% disp(strcat('MGRIT, iteration: ', num2str(MGcycle)));
% end
end
% store final error as well
if( sSol )
err(maxIT+1,:) = max( abs( u{1}(:,1:M(1):end)-uS(:,1:M(1):end) ), [], 1 );
end
% normalise error
% if( sSol )
% err = err/err(1);
% end
end
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
