package AMRElliptic;
import BoxTools.BoxLayout;
import BoxTools.BoxLayoutData;
import BoxTools.BoxedArray;
import BoxTools.Util;
import BoxTools.LayoutIterator;
import BoxTools.DataIndex;
import AMRTools.LevelFluxRegister;
import AMRTools.QuadCFInterp;
import java.lang.Math;
import ti.lang.Reduce;
/**
* AMRSolver provides a multi-level algorithm for elliptic problems. It should
* be implemented as a template class. The current implementation is the instance where the
* basic data type is double. A AMRSolver object is constructed based on a set
* of LevelOps, which defines the structure of the problem. Intermediate data
* structures and operators are reused to save resources.
* The details of the multigrid algorithms are available in Chombo Design Document.
* Usage: a LevelOp object must be declared local.
* @see BoxLayout,BoxLayoutData,BoxedArray,LayoutIterator,QuadCFInterp,CoarseAverage,ConstantFineInterp.
* @see Chombo Specification
* @version 1.1
* Modified on Jul 09, 2004. change of syntax of foreach, see BoxTools/foreachLoopTest.ti.
* on Aug 26, 2004. the operator BoxLayoutData.<<= constant covers ghost cells.
* @author Tong Wen, LBNL
* @since 1.0
*/
public class AMRSolver{
private static final int SPACE_DIM =Util.SPACE_DIM;
public int single numAMRVCycle=0;
LevelOp local single [] local m_LevelOps;
template BoxLayoutData > local single [] local phi, rhs, R, Rcopy, Rcoarse, e, phis, de;
template LevelFluxRegister local single [] local m_fluxRegisters;
private int single m_numLevels;
private double single m_epsilon;
Timer tmrCFInterp=new Timer();Timer tmrFluxReg=new Timer(); Timer tmrResidual=new Timer();Timer tmrMgRelax=new Timer();
//Timer tmrBSolver=new Timer();
public static Timer tmrUpdate=new Timer();public static Timer tmrAMRResidual=new Timer(); template QuadCFInterp cfi=new template QuadCFInterp();public static Timer tmrIncC=new Timer();public static Timer tmrIncf=new Timer();public static Timer tmrReflux=new Timer();public static Timer tmrGetflux=new Timer();public static Timer tmrAverage=new Timer();public static Timer tmrIpwc=new Timer();public static Timer tmrDeflux=new Timer(); public static Timer tmrAMRBLevel=new Timer();
public static Timer tmrCRRes=new Timer();public static Timer tmrCRCFInterp=new Timer();public static Timer tmrCRFlux=new Timer();
public AMRSolver( LevelOp local single [] local a_LevelOps){
m_LevelOps=a_LevelOps;
m_numLevels=a_LevelOps.length;
R=new template BoxLayoutData > local single [m_numLevels];
Rcopy=new template BoxLayoutData > local single [m_numLevels];
Rcoarse=new template BoxLayoutData > local single [m_numLevels-1];
e=new template BoxLayoutData > local single [m_numLevels];
phis=new template BoxLayoutData > local single [m_numLevels];
de=new template BoxLayoutData > local single [m_numLevels];
m_fluxRegisters= new template LevelFluxRegister local single [m_numLevels-1];
int single i;
BoxLayout local single boxLayout;
for (i=0;i >(boxLayout, 0);
Rcopy[i]=new template BoxLayoutData >(boxLayout, 0);
e[i]=new template BoxLayoutData >(boxLayout, 1);
phis[i]=new template BoxLayoutData >(boxLayout, 0);
de[i]=new template BoxLayoutData >(boxLayout, 1);
}
for (i=0;i >((m_LevelOps[i+1].m_layout).coarsen(m_LevelOps[i+1].m_refRatio), 0);
//System.out.println("AMRSolver:construct fluxRegisters at level "+i);
m_fluxRegisters[i]= new template LevelFluxRegister single (m_LevelOps[i].m_layout, m_LevelOps[i+1].m_layout,
m_LevelOps[i+1].m_refRatio ,m_LevelOps[i].m_Dx, m_LevelOps[i+1].m_domain);
//System.out.println("AMRSolver:construct fluxRegisters at level down!");
}
}
final local void computeAMRResidual(){
/*int single i;
i=m_numLevels-1;
computeAMRResidual(R[i],phi[i-1],phi[i],phi[i],rhs[i],i,true, true, true);
for (i=m_numLevels-2;i>0;i--){
computeAMRResidual(R[i],phi[i-1],phi[i],phi[i+1],rhs[i],i,true, true,true);
m_LevelOps[i+1].Average(R[i],R[i+1]);
}
computeAMRResidual(R[0],phi[0],phi[0],phi[1],rhs[0],0,true, true,true);
m_LevelOps[i+1].Average(R[0],R[1]);*/
int single i,j,k, dir, side;
LayoutIterator local LItr;
DataIndex local DI;
tmrCRCFInterp.start();
for (i=1;i0;i--) m_LevelOps[i].Average(R[i-1], R[i]); //correct order, but redundant for AMRVCycle
tmrCRFlux.stop();
}
//adding public for test purpose.
public final local single void computeAMRResidual(template BoxLayoutData > local single a_res,
template BoxLayoutData > local single a_phic,
template BoxLayoutData > local single a_phi,
template BoxLayoutData > local single a_phif,
template BoxLayoutData > local single a_rhs,
int single a_level, boolean single a_interpON,boolean single a_residualON,
boolean single a_refluxON){
tmrAMRResidual.start();
tmrCFInterp.start();
if (a_level>0 && a_interpON)
if (a_phi == a_phic)
m_LevelOps[a_level].CFInterp(a_phi);
else{
m_LevelOps[a_level].CFInterp(a_phi, a_phic);
}
tmrCFInterp.stop();
tmrResidual.start();
if (a_residualON)
m_LevelOps[a_level].residual(a_phi, a_rhs, m_LevelOps[a_level].m_Dx, a_res);
tmrResidual.stop();
//else
// m_fluxRegisters[a_level].deflux(a_res);
tmrFluxReg.start();
if (a_refluxON && a_level > local single [] local a_R,
int single a_normType){
LayoutIterator local LItr;
double [SPACE_DIM d] local TiArray;
//Point point; /* for change in the syntax of foreach */
double temp=0;
double single norm;
for (int single i=0;i > local single [] local a_R){
LayoutIterator local LItr;
double [SPACE_DIM d] local TiArray;
//Point point; /* for change in the syntax of foreach */
double max=0; double temp;
double single norm=0;
for (int single i=0;imax) max=temp;
}
LItr.advance();
}
//System.out.println("max norm at level "+i+" : "+max);
}
LItr=m_LevelOps[m_numLevels-1].m_layoutItr;
LItr.reset();
while(!LItr.isEnded()){
TiArray=a_R[m_numLevels-1][LItr.getDataIndex()].getLocalArray();
foreach (point in TiArray.domain()) {temp= Math.abs(TiArray[point]); if (temp>max) max=temp;}
LItr.advance();
}
//System.out.println("max norm at top level : "+m_LevelOps[m_numLevels-1].maxLocalNorm(a_R[m_numLevels-1])+" vs "+max);
norm=Reduce.max(max);
norm=broadcast norm from 0;
return norm;
}
public final local void solveAMR(template BoxLayoutData > local single [] local a_phi,
template BoxLayoutData > local single [] local a_rhs,
double single a_epsilon){
//(a_phi.length!=m_numLevels || a_rhs.length!=m_numLevels) cause segmentation fault.
if (a_epsilon<=0){
Util.printErrMsg("AMRSolver::solveAMR: Invalid inputs!");
return;
}
phi=a_phi; rhs=a_rhs;
Timer tmr1 =new Timer();Timer tmr2 =new Timer();Timer tmr3 =new Timer();
int single i,maxSteps;
//double single rhsNorm2=computeResidualNorm(rhs,2);
double single RNorm;
tmr1.start();computeAMRResidual();tmr1.stop();
//tmr2.start();RNorm=computeResidualNorm(R,2);tmr2.stop();
tmr2.start();RNorm=computeResidualMaxNorm(R);tmr2.stop();
if (numAMRVCycle==0)
maxSteps=Util.NUMAMRVCYCLE;
else
maxSteps=numAMRVCycle;
//maxSteps=1;
i=0;
Util.printMsg("Initial R norm is "+RNorm);
m_epsilon=a_epsilon*RNorm;
while (RNorm > m_epsilon && i0){
phis[l]<<=phi[l];
e[l]<<=0;
//tmrUpdate.start();e[l].exchange();m_LevelOps[l].applyHPhyBC(e[l]);tmrUpdate.stop();
tmrMgRelax.start();m_LevelOps[l].mgRelax(e[l], R[l]);tmrMgRelax.stop();
phi[l]+=e[l];
tmrUpdate.start();phi[l].exchange();m_LevelOps[l].applyPhyBC(phi[l]);tmrUpdate.stop();
//e[l-1]<<=0;e[l-1].exchange();m_LevelOps[l-1].applyPhyBC(e[l-1]);
tmrDeflux.start();m_fluxRegisters[l-1].deflux(R[l-1]);tmrDeflux.stop();
computeAMRResidual(R[l-1],phi[l-1],phi[l-1],phi[l],rhs[l-1],l-1,false,false,true);
computeAMRResidual(Rcopy[l],e[l],e[l],e[l],R[l],l,true, true,false);
tmrAverage.start();m_LevelOps[l].Average(R[l-1],Rcopy[l]);tmrAverage.stop();
AMRVCycleMG(l-1);
tmrIpwc.start();m_LevelOps[l].Ipwc(Rcopy[l], e[l-1]);tmrIpwc.stop();
e[l]+=Rcopy[l];
tmrUpdate.start();e[l].exchange();m_LevelOps[l].applyHPhyBC(e[l]);tmrUpdate.stop();
de[l]<<=0;
//tmrUpdate.start();de[l].exchange();m_LevelOps[l].applyHPhyBC(de[l]);tmrUpdate.stop();
computeAMRResidual(R[l],e[l-1],e[l],e[l],R[l],l,true, true,false);
tmrMgRelax.start();m_LevelOps[l].mgRelax(de[l], R[l]);tmrMgRelax.stop();
e[l]+=de[l];
//e[l].exchange();m_LevelOps[l].applyPhyBC(e[l]);//e[l-1] is used,but maybe redundant.
phi[l].setBySum(phis[l], e[l]);
tmrUpdate.start();phi[l].exchange();m_LevelOps[l].applyPhyBC(phi[l]);tmrUpdate.stop();
}
else{
tmrAMRBLevel.start();
e[0]<<=0;
//e[0].exchange();m_LevelOps[0].applyHPhyBC(e[0]);
/*tmrBSolver.start();m_LevelOps[0].bottomSmoother(e[0], R[0]);tmrBSolver.stop();*/
m_LevelOps[0].bottomSmoother(e[0], R[0]);
phi[0]+=e[0];
phi[0].exchange();m_LevelOps[0].applyPhyBC(phi[0]);
tmrAMRBLevel.stop();
}
//System.out.println("This is AMRSolver::AMRVCycleMG end at level "+a_level);
}
}