Overture and MetaChoas for CISM Developers

1. Overture and MetaChoas forCISM Developers M. Wiltberger With thanks to Bill Henshaw and Alan Sussman

2. What is Overture? • A Collection of C++ Classes that can be used to solve PDEs on overlapping grids • Key Features • High level interface for PDEs on adaptive and curvilinear grids • Provides a library of finite differences operators • Conservative/NonConservative • 2nd and 4Th order • Uses A++/P++ array class for serial parallel array operations • Extensive grid generation capablities

3. Solvers Oges, Ogmg, OverBlown Operators div, grad, bc's Grid Generator Ogen Adaptive Mesh Refinement Mappings (geometry) MappedGrid GridCollection MappedGridFunction GridCollectionFunction A++P++ array class Graphics (OpenGL) Data base (HDF) Boxlib (LBL) Overture: A toolkit for solving PDEs

4. C++ realArray u(10,10), v(10,10);Range I(8,8), J(8,8);v = 1;u(I,J) =0.25*(v(I+1,J)+v(I-1,J)+ v(I,J+1)+v(I,J-1)); Fortran do j=1,8 do i=1,8 u(i,j)=0.25*(v(i+1,j)+v(i-1,j) +v(i,j+1+v(i,j+1) enddoenddo A++ Array Class • Serial and Parallel Array operations • Easily passed to Fortran • Just recompile with different library for parallel code using MPI calls

5. extern "C" { void myfortran_( int & m, int & n, real & u); } realArray u(20,10); myfortran_( u.getLength(0), u.getLength(1), *u.getDataPointer() ); C++ subroutine myfortran( m,n,u ) real u(m,n) ... end Fortran Passing A++ arrays to Fortran

6. Each mapping has an optimized map and inverseMap map Mapping unit square domainDimension rangeDimension map( r,x,xr ) inverseMap( x,r,rx ) boundaryConditions singularities ... SquareMapping AnnulusMapping SphereMapping HyperbolicMapping EllipticTransform MatrixMapping RevolutionMapping etc. A Mapping defines a continuous transformation x r

7. Figure here MappedGrid gridIndexRange Mapping vertex vertexDerivative cellVolume faceNormal ... Mapping that defines the geometry grid points optionally computed geometry arrays jacobian derivatives A MappedGrid holds the grid points and other geometry info

8. base grids refinement grids GridCollection numberOfGrids operator [int g] refinementLevel[int l] ... access to a MappedGrid, g=0,1,.. GridCollection for a refinement level A GridCollection holds a list of MappedGrids

9. realMappedGridFunction numberOfComponents MappedGrid MappedGridOperators x,y,z,xx,... derivatives MappedGrid mg(mapping); realMappedGridFunction u(mg); u=1.; A MappedGridFunction holds solution values derived from an A++ realArray which implies it inherits all A++ operators scalar, vector, 2-tensor, .. Grid functions can be vertex-centered, cell-centered, face-centered etc.

10. realGridCollectionFunction realMappedGridFunction realMappedGridFunction GridCollection gc(...); Range all; realGridCollectionFunction u(gc,all,all,all,2); u=1.; for( int grid=0; grid<gc.numberOfGrids(); grid++ ) u[grid]=1.; A GridCollectionFunction is a list of MappedGridFunctions

11. The Operator classes define differential operators and boundary conditions. - 2nd/4th order approximations plus some conservative approximations MappedGrid mg(sphere); MappedGridOperators op(mg); floatMappedGridFunction u(mg), v(mg), w; v=u.x()+u.y(); w=u.laplacianCoefficients(); Compute the derivatives directly Form the sparse matrix for the differential operator Operators

12. GridCollection gc(...); realGridCollectionFunction u(gc), v(gc); GridCollectionOperators op(gc); v = u.x(); for( int grid=0; grid<gc.numberOfGrids(); grid++ ) v[grid]=u[grid].x(); for( grid=0; grid<gc.numberOfGrids(); grid++ ) xFortran_( *v[grid].getDataPointer(), *u[grid].getDataPointer(), ...); Operate at the collection level Operate at the single grid level Call Fortran Overture can be used at a number of levels

13. A sample Overture program Solve ut +aux + buy -(uxx+uyy) on overlaping grid int main() {CompositeGrid cg; // Create a Composite grid getFromADataBaseFile(cg,”mygrid.hdf”);floatCompositeGridFunction u(cg); u = 1.0; // Set initial conditionsCompositeGridOperators og(cg); // Create Operators u.setOperators(op); float t=0, dt=0.005, a=1., b=1., nu=0.1; for(int step=0;step<200;step++) {u+=dt*(-a*u.x()-b*u.y()+nu*(u.xx()+u.yy()); t+=dt;u.interpolate(); // Exchange info between gridsu.applyBoundaryConditions(dircichelt,allBoundaries); u.finishBoundaryConditions(); } return 0; }

14. Wave Equation Solver The Overture primer contains a short code showing how to write a solver for the wave equation. This example shows how a wave is able to pass through a grid interface with few artificial reflections.

15. Incompressible flow past two cylinders Fine grids are integrated implicitly, coarse grids explicitly.

16. What is MetaChoas? • A runtime meta-library that achieves direct data transfers between data structures managed by different parallel libraries • Runtime meta-library means that it interacts with data parallel libraries and languages used for separate applications (including MPI) • Can exchange data between separate (sequential or parallel) applications • Manage data transfers between different libraries in the same application • This often referred to as the MxN problem in parallel programming

17. How does MetaChaos work? • It starts with the Data Descriptor • Information about how the data is distributed across the processors • Usually supplied by the library/program developer • Sussman et al. are working on generalizing this process • MetaChaos then uses a linearization (LSA) of the data to be moved (the regions) to determine the optimal method to move data from set of regions in A (SA ) to a set of regions in B (SB) • Moving the data is a three step process LSA = llibX(SA)LSB = LSASB = l-1libY(LSB) • Only constraint on this operation is each region must have the same number of elements

18. A Simple Example: Wave Eq Using P++ #include <A++.h>main(int argc, char **argv) {Optimization_Manager::Initialize_Virtual_Machine("",iNPES,argc,argv); doubleArray daUnm1(iNumX+2,iNumY+2),daUn(iNumX+2,iNumY+2); doubleArray daUnp1(iNumX+2,iNumY+2); Index I(1,iNumX), J(1,iNumY); // Indices for computational domain for(j=1;j<iNumY+1;j++) { daUnm1(I,j) = sin(dW*dTime + (daX(I)*2*dPi)/dLenX); daUn(If,j) = sin(dW*0 + (daX(If)*2*dPi)/dLenX); } // Apply BC Omitteed for space // Evole a step forward in time for(i=1;i<iNSteps;i++) { daUnp1(I,J) = ((dC*dC*dDT*dDT)/(dDX*dDX))* (daUn(I-1,J)-2*daUn(I,J)+daUn(I+1,J)) + 2*daUn(I,J) - daUnm1(I,J); // Apply BC Omitted for space }Optimization_Manager::Exit_Virtual_Machine();}

19. Split into two using MetaChaos #include <A++.h>main(int argc, char **argv) {Optimization_Manager::Initialize_Virtual_Machine("",NPES,argc,argv);this_pgme = InitPgme(pgme_name,NPES); other_pgme = WaitPgme(other_pgme_name,NPES_other); Sync2Pgme(this_pgme,other_pgme); BP_set = Alloc_setOfRegion(); left[0] = 4; right[0] = 4; stride[0] = 1; left[1] = 5; right[1] = 5; stride[0] = 1; reg = Alloc_R_HDF(DIM,left,right,stride); Add_Region_setOfRegion(reg,BP_Set); BP_da = getPartiDescriptor(&daUn); sched = ComputeScheduleForSender(…,BP_da,BP_set,…); for(i=1;i<iNSteps;i++) { daUnp1(I,J) = ((dC*dC*dDT*dDT)/(dDX*dDX))* (daUn(I-1,J)-2*daUn(I,J)+daUn(I+1,J)) + 2*daUn(I,J) - daUnm1(I,J);iDataMoveSend(other_pgme,sched,daUn,getLocalArray().getDataPointer); iDataMoveRecv(other_pgme,sched,daUn,getLocalArray().getDataPointer); Sync2Pgme(this_pgme,other_pgme); }Optimization_Manager::Exit_Virtual_Machine();}

20. How CISM will use Overture and MetaChoas