[overture] Re: deforming grid PDE solver
- From: jacob@xxxxxxxxxxxx
- To: overture@xxxxxxxxxxxxx
- Date: Thu, 26 Jul 2012 13:11:22 -0400
Hi Bill,
Here is a first stab: I just took deform.C and tried to tack on a
diffusion equation solver on top. The current code gives rather
spectacularly bad results, which is a little confusing to me, as my much
more complicated fluid solver, which to my understanding works more or
less the same, does not fail nearly this badly...Here it is:
#include "Ogen.h"
#include "PlotStuff.h"
#include "NurbsMapping.h"
#include "HyperbolicMapping.h"
#include "Ogshow.h"
#include "CompositeGridOperators.h"
//
// Deforming Grid Example:
// o
// o interpolate a grid function, update the interpolant for the new grid
// o save solutions in a show file
//
int
main(int argc, char *argv[])
{
Overture::start(argc,argv); // initialize Overture
Mapping::debug=0;
int numberOfSteps=50; // take small time-steps for the grid
int plotOption=true;
if( argc > 1 )
{ // look at arguments for "noplot"
aString line;
int len=0;
for( int i=1; i<argc; i++ )
{
line=argv[i];
if( line=="noplot" )
plotOption=false;
else if( len=line.matches("-numSteps=") )
{
sScanF(line(len,line.length()-1),"%i",&numberOfSteps);
printF("Setting numberOfSteps=%i\n",numberOfSteps);
}
}
}
aString nameOfOGFile="iceCircle.hdf"; // use this grid to start
// cout << "Enter the name of the (old) overlapping grid file:" << endl;
// cin >> nameOfOGFile;
// Create two CompositeGrid objects, cg[0] and cg[1]
CompositeGrid cg[2];
getFromADataBase(cg[0],nameOfOGFile); // read cg[0] from a
data-base file
cg[1]=cg[0]; // copy cg[0] into cg[1]
// Deform component grid 2 (do this by changing the mapping)
int gridToDeform=2;
// The grid we deform will be a Hyperbolic Mapping, we need 2 of them:
HyperbolicMapping transform[2];
NurbsMapping *startCurve = new NurbsMapping [2]; // start curve
representing the ice surface
const int n0=13;
realArray x0(n0,2), x1(n0,2), x2(n0,2);
real rad=.5, theta;
// x0: points on a arc
x0(0,0)=0.; x0(0,1)=1.*rad;
theta=Pi* 5./180.; x0(1,0)=-rad*sin(theta); x0(1,1)=rad*cos(theta);
theta=Pi*10./180.; x0(2,0)=-rad*sin(theta); x0(2,1)=rad*cos(theta);
theta=Pi*15./180.; x0(3,0)=-rad*sin(theta); x0(3,1)=rad*cos(theta);
theta=Pi*20./180.; x0(4,0)=-rad*sin(theta); x0(4,1)=rad*cos(theta);
theta=Pi*50./180.; x0(5,0)=-rad*sin(theta); x0(5,1)=rad*cos(theta);
theta=Pi*90./180.; x0(6,0)=-rad*sin(theta); x0(6,1)=rad*cos(theta);
theta=Pi*130./180.;x0(7,0)=-rad*sin(theta); x0(7,1)=rad*cos(theta);
x0( 8,0)=x0(4,0); x0( 8,1)=-x0(4,1);
x0( 9,0)=x0(3,0); x0( 9,1)=-x0(3,1);
x0(10,0)=x0(2,0); x0(10,1)=-x0(2,1);
x0(11,0)=x0(1,0); x0(11,1)=-x0(1,1);
x0(12,0)=x0(0,0); x0(12,1)=-x0(0,1);
// a second start curve (deformed version)
x1=x0;
x1(5,0)=-.7; x1(5,1)= .5;
x1(6,0)=-.6; x1(6,1)= 0.;
x1(7,0)=-.7; x1(7,1)= -.5;
startCurve[0].interpolate(x0);
startCurve[1].interpolate(x0);
for( int i=0; i<=1; i++ )
{
const bool isSurfaceGrid=false, init=true;
transform[i].setSurface(startCurve[i],isSurfaceGrid,init);
transform[i].setGridDimensions(axis1,61);
transform[i].setParameters(HyperbolicMapping::distanceToMarch,.25);
transform[i].setParameters(HyperbolicMapping::linesInTheNormalDirection,9);
transform[i].setShare(0,1,1);
transform[i].generate();
}
// Replace the mapping of the component grid that we want to move:
for( int i=0; i<=1; i++ )
{
cg[i][gridToDeform].reference(transform[i]);
cg[i].update();
}
// now we destroy all the data on the new grid -- it will be shared with
the old grid
// this is not necessary to do but it will save space
cg[1].destroy(CompositeGrid::EVERYTHING);
// we tell the grid generator which grids have changed
LogicalArray hasMoved(cg[0].numberOfGrids());
hasMoved = LogicalFalse;
hasMoved(gridToDeform) = LogicalTrue; // Only this grid will change
char buff[80];
//set boundary conditions in the proper locations
const int DIRICHLET = 1;
cg[0][0].boundaryCondition()(0,0) = DIRICHLET;
cg[0][0].boundaryCondition()(1,0) = DIRICHLET;
cg[0][0].boundaryCondition()(0,1) = DIRICHLET;
cg[0][0].boundaryCondition()(1,1) = DIRICHLET;
cg[0][1].boundaryCondition()(0,0) = 0;
cg[0][1].boundaryCondition()(1,0) = 0;
cg[0][1].boundaryCondition()(0,1) = DIRICHLET;
cg[0][1].boundaryCondition()(1,1) = 0;
cg[0][2].boundaryCondition()(0,0) = 0;
cg[0][2].boundaryCondition()(1,0) = 0;
cg[0][2].boundaryCondition()(0,1) = DIRICHLET;
cg[0][2].boundaryCondition()(1,1) = 0;
PlotStuff ps(plotOption,"Deforming Grid Example"); // for plotting
PlotStuffParameters psp;
// Here is the overlapping grid generator
Ogen gridGenerator(ps);
// update the initial grid, since the above reference destroys the mask
gridGenerator.updateOverlap(cg[0]);
// Here is an interpolant
Interpolant & interpolant = *new Interpolant(cg[0]);
interpolant.incrementReferenceCount(); //why do we have to do this here?
CompositeGridOperators operators(cg[0]); //
operators for the CompositeGrid
realCompositeGridFunction u(cg[0]), F(cg[0]);
u = 0.;
F = 1.;
u.setOperators(operators);
F.setOperators(operators);
//diffusion term
const real nu = 1;
//we'll solve the diffusion equation
/******************
u_t = nu*(u_{xx} + u_{yy}) + 1
u(0) = 0
******************/
//with zero dirichlet boundary conditions on all boundaries
//the first version will be forward Euler explicit time stepping
// Here is a show file
Ogshow show("deform.show");
show.setIsMovingGridProblem(true);
//the time step size for the grid deformation
real delta=1./numberOfSteps; //moves the time-step
int nx = 81; //the grid dimensions for CFL
//the time-step size for the forward Euler solves
const int numSubSteps = std::ceil(3*nx*nx);
const real dt = 1./numSubSteps; //viscous time-step restriction for the
forward Euler solves
std::cout << "numSubSteps: " << numSubSteps << std::endl;
// ---- Deform the grid in a periodic fashion.----
for (int i=1; i<=numberOfSteps; i++) {
int newCG = i % 2; // new grid
int oldCG = (i+1) % 2; // old grid
ps.erase();
psp.set(GI_TOP_LABEL,sPrintF(buff,"Solution at step=%i",i)); // set
title
//PlotIt::plot(ps,cg[oldCG],psp); // plot the current
overlapping grid
PlotIt::contour(ps,u,psp);
if( i>10 )
psp.set(GI_PLOT_THE_OBJECT_AND_EXIT,TRUE); // set this to
run in "movie" mode (after first plot)
ps.redraw(TRUE);
//real deltaOmega=1./20.;
real deltaOmega = delta;
real omega=pow(sin(Pi*i*deltaOmega),2.); // omega varies in the
interval [0,1]
x2 = (1.-omega)*x0 + omega*x1;
startCurve[newCG].interpolate(x2); // form the new start curve
(NurbsMapping)
const bool isSurfaceGrid=false;
const bool init=false; // this means keep existing hype parameters
such as distanceToMarch, linesToMarch etc.
transform[newCG].setSurface(startCurve[newCG],isSurfaceGrid,init); //
supply a new start curve
// ** generate the new grid with the hyperbolic grid generator **
transform[newCG].generate();
// Update the overlapping grid newCG, starting with and sharing data
with oldCG.
Ogen::MovingGridOption option = Ogen::useOptimalAlgorithm;
int useFullAlgorithmInterval=1; // 10;
if( i% useFullAlgorithmInterval == useFullAlgorithmInterval-1 )
{
cout << "\n +++++++++++ use full algorithm in updateOverlap
+++++++++++++++ \n";
option=Ogen::useFullAlgorithm;
}
gridGenerator.updateOverlap(cg[newCG], cg[oldCG], hasMoved, option );
cg[newCG][0].boundaryCondition()(0,0) = DIRICHLET;
cg[newCG][0].boundaryCondition()(1,0) = DIRICHLET;
cg[newCG][0].boundaryCondition()(0,1) = DIRICHLET;
cg[newCG][0].boundaryCondition()(1,1) = DIRICHLET;
cg[newCG][1].boundaryCondition()(0,0) = 0;
cg[newCG][1].boundaryCondition()(1,0) = 0;
cg[newCG][1].boundaryCondition()(0,1) = DIRICHLET;
cg[newCG][1].boundaryCondition()(1,1) = 0;
cg[newCG][2].boundaryCondition()(0,0) = 0;
cg[newCG][2].boundaryCondition()(1,0) = 0;
cg[newCG][2].boundaryCondition()(0,1) = DIRICHLET;
cg[newCG][2].boundaryCondition()(1,1) = 0;
cg[newCG].update();
//move the solution and operators to the new grid [!?]
interpolant.updateToMatchGrid(cg[newCG]);
operators.updateToMatchGrid(cg[newCG]);
u.updateToMatchGrid(cg[newCG]);
interpolant.interpolate(u);
F.updateToMatchGrid(cg[newCG]);
interpolant.interpolate(F);
//clamp u to zero on the boundaries
for(int grid = 0; grid <= 2; grid++) { //ok I know there are exactly 3
grids
//MappedGrid& mg = cg[newCG][grid];
for(int axis = 0; axis <= 1; axis++) {
for(int side = 0; side <= 1; side++) {
if( cg[newCG][grid].boundaryCondition()(side,axis) ==
DIRICHLET) {
Index I1,I2,I3;
getBoundaryIndex(cg[newCG][grid].gridIndexRange(),side,axis,I1,I2,I3);
u[grid](I1,I2,I3) = 0.;
}
}
}
}
u.applyBoundaryCondition(0,BCTypes::dirichlet,DIRICHLET,0.);
u.finishBoundaryConditions();
//take the requisite number of forward Euler time-steps to catch up to
the grid
for(int k = 0; k < numSubSteps; k++) {
u += dt*( nu*(u.xx() + u.yy()) + F);
}
interpolant.interpolate(u);
// save the result in a show file, every fourth step
if( (i % 4) == 1 )
{
show.startFrame();
show.saveComment(0,sPrintF(buff,"Solution at step = %i",i));
show.saveSolution(u);
}
}
printF("Results saved in deform.show, use Overture/bin/plotStuff to view
this file\n");
printF("deform:
interpolant.getReferenceCount=%i\n",interpolant.getReferenceCount());
if( interpolant.decrementReferenceCount()==0 )
{
printF("deform: delete Interpolant\n");
delete &interpolant;
}
Overture::finish();
return 0;
}
Best,
Jacob
On Thu, July 26, 2012 9:02 am, Bill Henshaw wrote:
> Hi Jacob,
> My thought was to develop a relatively simple example code
> that we could add to the Overture distribution for other
> people to use. If you want to make a first attempt then
> I can make adjustments.
>
> ....Bill
>
>
> On 07/25/2012 09:49 PM, jacob@xxxxxxxxxxxx wrote:
>> Hi Bill,
>>
>> Hmm I agree, though I had thought move2.C was a bit complicated and
>> wasn't
>> sure what exactly it was doing, so I more used deform.C to make my guess
>> at how to deform the grids...But yeah I can take a pass at making a test
>> code which is a bit simpler...maybe I'll try advection-diffusion first
>> with explicit and then with implicit stepping, since updating the solver
>> should "probably" be trickier. For the deforming boundary I'll probably
>> first use some kind of spline that deforms in a pre-defined way and try
>> to
>> build a HyperbolicMapping off of it...if that doesn't work I'll
>> back-pedal
>> and try something simpler...
>>
>> Best,
>> Jacob
>>
>> On Wed, July 25, 2012 11:47 pm, Bill Henshaw wrote:
>>> Hi Jacob,
>>> Thanks for the code segments. You have made some good progress,
>>> although
>>> I do see some problems with what you are doing.
>>>
>>> Recall:
>>> GOLDEN RULE 3 of programming -- "It always saves time in
>>> the long run to have a simple test code".
>>>
>>> Corollary of Golden Rule 3: you may hate this rule but unfortunately
>>> it is really true.
>>>
>>> Unfortunately I have not created a simple deforming grid
>>> PDE solver code. I think that the first step is for us to create
>>> such a code. I think the easiest way to do this is to merge
>>> deform.C and move2.C so that we solve an advection diffusion
>>> equation with twilight-zone on a deforming grid. Perhaps you
>>> want to make a first attempt at doing this? Otherwise I can
>>> try to put something together but I'm not sure when I can get
>>> to it.
>>>
>>>
>>> ...Bill
>>>
>>>
>>>
>>
>>
>>
>
>
Other related posts: