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
.
