[overture] Re: interpolation error?
- From: jacob@xxxxxxxxxxxx
- To: overture@xxxxxxxxxxxxx
- Date: Wed, 25 Jul 2012 14:11:15 -0400
Hi Bill,
Well the grids are generated in my C++ code, but here what I guess should
be the relevant code snippets. The test that is breaking doesn't move the
grid, but it updates the CompositeGrid as if I did, so it might have
something to do with that update. Here is the code:
//somewhere up here the grid is first made:
Mapping* domain;
domain = new AnnulusMapping(.4,R);
domain->setGridDimensions(axis1,std::floor(Ntheta) );
domain->setGridDimensions(axis2,std::floor(Nr*.5)); //
domain->setBoundaryCondition(0,1,0);
Mapping& Domain = *domain;
//this is a class which stores the interface in a format which is useful
for the problem at hand
//the constructor defaults to initialize it into a circle of radius one.
ThetaLMapping curve(50);
curve.setGridDimensions(0,Ntheta);
curve.setIsPeriodic(0,Mapping::functionPeriodic);
Range all;
//I found that the only reliable way to get the code working is to
translate my object into a spline
SplineMapping spline(2);
spline.setGridDimensions(0,Ntheta);
spline.setIsPeriodic(0,Mapping::functionPeriodic);
realArray r(Ntheta,1,1);
realArray x(Ntheta,2,1);
for(int k = 0; k <=r.getBound(0); k++) {
r(k) = ((double)k)/(r.getBound(0)+1);
x(k,0) = 0;
x(k,1) = 0;
}
curve.map(r,x);
spline.setPoints(x(all,0),x(all,1));
spline.setBoundaryCondition( 0,0,-1 );
spline.setBoundaryCondition( 1,0,-1 );
//build a body-fit grid around the spline
HyperbolicMapping hyper(spline);
hyper.setIsPeriodic(0,Mapping::functionPeriodic);
hyper.setParameters(HyperbolicMapping::distanceToMarch,BLthickness);
hyper.setParameters(HyperbolicMapping::linesInTheNormalDirection,NBL,HyperbolicMapping::forwardDirection);
hyper.generate();
hyper.setBoundaryCondition(1,1,0);
PlotStuff gi(true,"hyperbolic map");
MappingInformation mapInfo;
mapInfo.graphXInterface=&gi;
mapInfo.mappingList.addElement(Domain);
mapInfo.mappingList.addElement(hyper);
const int numberOfDimensions=2, numberOfGrids=2;
IntegerArray mapList(numberOfGrids);
mapList(0)=0;
mapList(1)=1;
//makes the composite grid
CompositeGrid compGrid;
Ogen ogen(gi);
ogen.buildACompositeGrid(compGrid,mapInfo,mapList);
ogen.updateOverlap(compGrid);
compGrid.update();
const int ABC=1, NOSLIP=2,PERIODIC=-1;
//zero should stand for extrapolation BCs
compGrid[0].boundaryCondition()(0,1)=0;
compGrid[0].boundaryCondition()(1,1)=ABC; //artificial boundary condition
for outflow/inflow
compGrid[1].boundaryCondition()(0,1)=NOSLIP;
compGrid[1].boundaryCondition()(1,1)=0;
compGrid[1].boundaryCondition()(0,0)=PERIODIC;
compGrid[1].boundaryCondition()(1,0)=PERIODIC;
realCompositeGridFunction force(compGrid), psi(compGrid), omega(compGrid),
u(compGrid), v(compGrid), f(compGrid); //the latter holding the body
forces and various RHS for solves.
int stencilSize=int( pow(3,compGrid.numberOfDimensions())+1 ); // add 1
for interpolation equations wtf?
realCompositeGridFunction coeffLap(compGrid,stencilSize,all,all,all);
coeffLap.setIsACoefficientMatrix(TRUE,stencilSize,1);
coeffLap = 0.;
CompositeGridOperators operators(compGrid);
operators.setStencilSize(stencilSize);
operators.setOrderOfAccuracy(2);
omega.setOperators(operators);
u.setOperators(operators);
v.setOperators(operators);
f.setOperators(operators);
force.setOperators(operators);
psi.setOperators(operators);
omega = 0.;
force = 0.;
f = u = v = psi = 0.;
Interpolant & interpolant = *new Interpolant(compGrid); //
do this instead for now.
interpolant.setImplicitInterpolationMethod(Interpolant::iterateToInterpolate);
coeffLap.setOperators(operators);
coeffLap = operators.laplacianCoefficients();
/*************
the artificial BCs are also set here etc, but let's skip that
**************/
//the Oges object is defined like this:
Oges PoissonSolver(compGrid); // create a solver
PoissonSolver.setCoefficientArray( coeffLap ); // supply coefficients
/****************************************
Then later in the time stepping:
*****************************************/
//here is where I would evolve the curve object, but I don't. For this
test //it is just stays a circle, but most of the other code doesn't know
that
//update the grids and composite grid functions
//2. First we have to re-build the hyperbolic mapping and composite grid
CompositeGrid oldGrid = compGrid;
spline.setIsPeriodic(0,Mapping::functionPeriodic);
spline.setGridDimensions(0,Ntheta);
realArray r2(shearStress.size(),1,1);
realArray x2(shearStress.size(),2,1);
for(int k = 0; k <=r2.getBound(0); k++) {
r2(k) = ((double)k)/(r2.getBound(0)+1);
x2(k,0) = 0;
x2(k,1) = 0;
}
curve.map(r2,x2);
spline.setPoints(x2(all,0),x2(all,1));
spline.setBoundaryCondition( 0,0,-1 );
spline.setBoundaryCondition( 1,0,-1 );
HyperbolicMapping hyperReplace(spline);
hyperReplace.setIsPeriodic(0,Mapping::functionPeriodic);
hyperReplace.setParameters(HyperbolicMapping::distanceToMarch,BLthickness);
hyperReplace.setParameters(HyperbolicMapping::linesInTheNormalDirection,NBL,HyperbolicMapping::forwardDirection);
hyperReplace.generate();
hyperReplace.setBoundaryCondition(1,1,0);
hyper=hyperReplace;
std::cout << "rebuilt hyperbolic mapping" << std::endl;
compGrid.update();
compGrid[0].boundaryCondition()(0,1)=0;
compGrid[0].boundaryCondition()(0,0)=PERIODIC;
compGrid[0].boundaryCondition()(1,0)=PERIODIC;
compGrid[0].boundaryCondition()(1,1)=ABC;
compGrid[1].boundaryCondition()(0,1)=NOSLIP;
compGrid[1].boundaryCondition()(1,1)=0;
compGrid[1].boundaryCondition()(0,0)=PERIODIC;
compGrid[1].boundaryCondition()(1,0)=PERIODIC;
MappingInformation mapInfo2;
mapInfo2.graphXInterface=&gi;
mapInfo2.mappingList.addElement(Domain);
mapInfo2.mappingList.addElement(hyper);
IntegerArray mapList2(numberOfGrids);
mapList2(0)=0;
mapList2(1)=1;
LogicalArray hasMoved(2);
hasMoved = LogicalTrue;
ogen.updateOverlap(compGrid,oldGrid,hasMoved,Ogen::useFullAlgorithm);
compGrid.update();
std::cout << "rebuilt composite grid" << std::endl;
//3. next we have to interpolate the previous time-step to the new grid
operators.updateToMatchGrid(compGrid);
//move all of the data to the new grids?
omega.setOperators(operators);
u.setOperators(operators);
v.setOperators(operators);
f.setOperators(operators);
force.setOperators(operators);
psi.setOperators(operators);
interpolant.updateToMatchGrid(compGrid);
force.updateToMatchGrid(compGrid);
force.interpolate();
psi.updateToMatchGrid(compGrid);
psi.interpolate();
omega.updateToMatchGrid(compGrid);
omega.interpolate();
u.updateToMatchGrid(compGrid);
u.interpolate();
v.updateToMatchGrid(compGrid);
v.interpolate();
f.updateToMatchGrid(compGrid);
f.interpolate();
//4. update the matrix for the linear solver...I'm not sure but it looks
like I have to reset Oges
coeffLap.updateToMatchGrid(compGrid);
coeffLap.setOperators(operators);
coeffLap.setIsACoefficientMatrix(TRUE,stencilSize,1);
coeffLap = 0.;
coeffLap = operators.laplacianCoefficients();
coeffLap.applyBoundaryConditionCoefficients(0,0,BCTypes::dirichlet, NOSLIP);
coeffLap.applyBoundaryConditionCoefficients(0,0,BCTypes::mixed,ABC,bcParams);
//set mixed conditions.
coeffLap.applyBoundaryConditionCoefficients(0,0,BCTypes::extrapolate,NOSLIP);
// extrap ghost line
coeffLap.applyBoundaryConditionCoefficients(0,0,BCTypes::extrapolate,0);
// extrap ghost line
coeffLap.finishBoundaryConditions();
PoissonSolver.setGrid(compGrid);
PoissonSolver.updateToMatchGrid(compGrid);
PoissonSolver.setCoefficientArray(coeffLap);
Thanks!
Best,
Jacob
On Wed, July 25, 2012 11:34 am, Bill Henshaw wrote:
> Hi Jacob,
> Can you send me the command file to generate the grid?
>
> ...Bill
>
> On 07/25/2012 08:20 AM, jacob@xxxxxxxxxxxx wrote:
>> Hi,
>>
>> I'm working on a composite grid which is basically two overlapping
>> annular
>> regions and at high resolution it breaks down with a lot of errors like
>> this:
>>
>> Interpolant:WARNING: no convergence in implicit interpolation iteration,
>> 25 its, max residual=5.136851e+207
>>
>> The eccentricity of the elements doesn't seem that bad, especially where
>> the overlap is, so I'm not sure what could be causing this...Thanks!
>>
>> Best,
>> Jacob
>>
>>
>>
>>
>
>
Other related posts:
