vtf-logo

AMRELCGFMSolver.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2003-2007 California Institute of Technology
00004 // Ralf Deiterding, ralf@amroc.net
00005 
00006 #ifndef AMROC_ELC_GFMSOLVER_H
00007 #define AMROC_ELC_GFMSOLVER_H
00008 
00016 #include "AMRCoupledGFMSolver.h"
00017 #include "AMRGFMInterpolation.h"
00018 #include "CPTLevelSet.h"
00019 #include "Timing.h"
00020 #include "elc/EulerianCommBoundary.h"
00021 #include "elc/EulerianCommShell.h"
00022 #include <vector>
00023 
00024 
00031 template <class VectorType, class FixupType, class FlagType, int dim>
00032 class AMRELCGFMSolver : public AMRCoupledGFMSolver<VectorType,FixupType,FlagType,dim> {
00033   typedef AMRCoupledGFMSolver<VectorType,FixupType,FlagType,dim> base;
00034   typedef typename VectorType::InternalDataType DataType;
00035 public:
00036   typedef typename base::integrator_type integrator_type;
00037   typedef typename base::initial_condition_type initial_condition_type;
00038   typedef typename base::boundary_conditions_type boundary_conditions_type;
00039   typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
00040 
00041   typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,dim> interpolation_type;
00042   typedef typename interpolation_type::point_type point_type;
00043   typedef typename elc::EulerianComm<dim,DataType> eul_comm_type;
00044   typedef typename elc::EulerianCommShell<dim,DataType>    eul_comm_pressure_on_face_type;
00045   typedef typename elc::EulerianCommBoundary<dim,DataType> eul_comm_pressure_on_node_type;
00046   typedef typename geom::BBox<dim,DataType> bbox_type;
00047   typedef CPTLevelSet<DataType,dim> cpt_type;
00048 
00049   AMRELCGFMSolver(integrator_type& integ, 
00050                   initial_condition_type& init,
00051                   boundary_conditions_type& bc) : base(integ, init, bc), _CoupleGFM(-1),
00052     _SendDataLocation(0), _ThinStructure(0), _LengthConversionFactor(1.), _PressureConversionFactor(1.) {      
00053     _Interpolation = new interpolation_type(*this);
00054     _eulComm = (eul_comm_type*) 0; 
00055   }
00056 
00057   virtual ~AMRELCGFMSolver() {
00058     delete _Interpolation;
00059     if (_eulComm) delete _eulComm;
00060   }
00061     
00062   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00063     base::register_at(Ctrl,prefix);
00064     RegisterAt(base::LocCtrl,"DataLocation",_SendDataLocation);
00065     RegisterAt(base::LocCtrl,"ThinStructure",_ThinStructure);
00066     RegisterAt(base::LocCtrl,"LengthFactor",_LengthConversionFactor);
00067     RegisterAt(base::LocCtrl,"PressureFactor",_PressureConversionFactor);
00068   } 
00069   virtual void register_at(ControlDevice& Ctrl) {
00070     base::register_at(Ctrl);
00071   }
00072 
00073   virtual void finish() {
00074     if (_eulComm) {
00075       delete _eulComm;
00076       _eulComm = (eul_comm_type*) 0;
00077     }
00078     base::finish();
00079   }
00080 
00081   virtual void SetupData() {
00082     base::SetupData();
00083     _Interpolation->SetupData(base::PGH(), base::NGhosts());
00084     if (_SendDataLocation!=0) _SendDataLocation=1;
00085   }
00086 
00087   void SetupInterComm(const int solid_nodes, const int solid_root) {
00088     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00089     if (SendDataLocation()==0)
00090       _eulComm = new eul_comm_pressure_on_face_type(MPI_COMM_WORLD, base::Comm(), solid_nodes,
00091                                                     solid_root, elc::GlobalIdentifiers);
00092     else
00093       _eulComm = new eul_comm_pressure_on_node_type(MPI_COMM_WORLD, base::Comm(), solid_nodes,
00094                                                     solid_root, elc::GlobalIdentifiers);
00095 #ifdef DEBUG_PRINT_ELC
00096     ( comm_service::log() << "*** EulerianComm created: SolidNodes=" << solid_nodes 
00097       << " SolidRoot=" << solid_root << "\n" ).flush();
00098 #endif
00099   }
00100 
00101   virtual void ReceiveBoundaryData(bool FullDomain=false) {
00102     START_WATCH
00103       BBox wholebb;
00104 
00105       if (!FullDomain) {
00106         int Level = 0;
00107         int Time = CurrentTime(base::GH(),Level);
00108         forall (base::U(),Time,Level,c)      
00109           wholebb += base::U().boundingbbox(Time,Level,c);
00110         end_forall
00111       }
00112       else 
00113         wholebb = base::GH().wholebbox();
00114 
00115       DCoords lbcorner = base::GH().worldCoords(wholebb.lower(), wholebb.stepsize());
00116       DCoords ubcorner = base::GH().worldCoords(wholebb.upper()+wholebb.stepsize(), 
00117                                                 wholebb.stepsize());
00118       if (LengthConversionFactor()!=1.) {
00119         lbcorner /= LengthConversionFactor();
00120         ubcorner /= LengthConversionFactor();
00121       }
00122       bbox_type lregion(lbcorner(), ubcorner());
00123             
00124 #ifdef DEBUG_PRINT_ELC
00125       ( comm_service::log() << "*** ELC::ReceiveBoundary, " << lregion ).flush();
00126 #endif
00127 
00128       START_WATCH
00129         _eulComm->receiveMesh(lregion);
00130         _eulComm->waitForMesh();
00131       END_WATCH(ELC_RECEIVEBOUNDARY)
00132 
00133       ModifyELCReceiveData(_eulComm);
00134 
00135 #ifdef DEBUG_PRINT_ELC
00136       ( comm_service::log() << " done.\n" ).flush();
00137 #endif
00138       int Npoints = _eulComm->getNumberOfNodes();
00139       point_type* positions = (point_type *)(_eulComm->getPositionsData());
00140       point_type* velocities = (point_type *)(_eulComm->getVelocitiesData());
00141       register int n;
00142       if (LengthConversionFactor()!=1.) {
00143         for (n=0; n<Npoints; n++) {
00144           positions[n] *= static_cast<DataType>(LengthConversionFactor());
00145           velocities[n] *= static_cast<DataType>(LengthConversionFactor());
00146         }
00147       }
00148 #ifdef DEBUG_PRINT_ELC_VALUES
00149       ( comm_service::log() << "*** Positions received from ELC ***\n" ).flush();
00150       for (n=0; n<Npoints; n++) 
00151         ( comm_service::log() << positions[n] << "\n" ).flush();    
00152 #endif
00153     END_WATCH(FLUID_CPL_RECEIVE_OVERHEAD)
00154 
00155       (reinterpret_cast<cpt_type* >(base::GFM(CoupleGFM()).GetLevelSetP()))->SetBrep(_eulComm->getNumberOfNodes(),
00156                  _eulComm->getPositionsData(), _eulComm->getNumberOfFaces(),
00157                  _eulComm->getConnectivitiesData());
00158   }
00159 
00160   virtual void SendBoundaryData() {
00161 #ifdef DEBUG_PRINT
00162     ( comm_service::log() << "*** SendBoundaryData()\n" ).flush();
00163 #endif
00164     START_WATCH
00165 #ifdef DEBUG_PRINT
00166       ( comm_service::log() << "*** Geometry construction()\n" ).flush();
00167 #endif
00168       DataType distance = -base::GFM(CoupleGFM()).Boundary().Cutoff();
00169       
00170       int Time = CurrentTime(base::GH(),base::CouplingLevel);
00171       DataType dxmin = 1.e37;
00172       forall(base::U(),Time,base::CouplingLevel,c)
00173         DCoords dx = base::GH().worldStep(base::U()(Time,base::CouplingLevel,c).stepsize());
00174         for (register int d=0; d<base::Dim(); d++) 
00175           if (dx(d)<dxmin) dxmin=dx(d);
00176       end_forall 
00177       if (distance>dxmin*base::NGhosts()) {
00178         int ghn = static_cast<int>(distance/dxmin + 1);
00179         std::cerr << "\n\nToo few ghostcells for Cutoff=" << -distance 
00180                   << ". Set GhostCells>=" << ghn << " or set Cutoff<=" 
00181                   << dxmin*base::NGhosts() << ".\n";
00182         std::cerr << "Aborting programm...\n" << std::flush;
00183         std::exit(-1);
00184       }
00185 
00186       _eulComm->initializePressure();
00187 
00188       register int n;
00189       int Npoints;
00190 
00191       double *press_data = _eulComm->getPressuresData();
00192       const int* con = _eulComm->getConnectivitiesData();
00193       const point_type *pos = reinterpret_cast<const point_type *>(_eulComm->getPositionsData());
00194       if (_ThinStructure>0 || distance>1.e-12)
00195         _eulComm->computeFaceNormals();
00196       if (_ThinStructure>0 || SendDataLocation()!=1 || distance>1.e-12)
00197         _eulComm->computeFaceCentroids();
00198       const point_type *norm = reinterpret_cast<const point_type *>(_eulComm->getFaceNormalsData());
00199       const point_type *cent = reinterpret_cast<const point_type *>(_eulComm->getFaceCentroidsData());
00200 
00201       if (SendDataLocation()==1)
00202         Npoints = _eulComm->getNumberOfNodes();
00203       else
00204         Npoints = _eulComm->getNumberOfFaces();
00205 
00206       if (_ThinStructure>0) {
00207         point_type* positions = new point_type[2*Npoints];
00208         point_type* centroids = new point_type[2*Npoints];
00209         if (SendDataLocation()==1) {
00210           DataType* face_count = new DataType[Npoints];
00211           point_type* point_normal = new point_type[Npoints];
00212           START_WATCH
00213             for (n=0; n<Npoints; n++) {
00214               face_count[n] = static_cast<DataType>(0.);
00215               point_normal[n] = static_cast<point_type>(0.);
00216             }
00217               
00218             int Nfaces = _eulComm->getNumberOfFaces();
00219             for (n=0; n<Nfaces; n++) 
00220               for (register int d=0; d<dim; d++) {
00221                 point_normal[con[dim*n+d]] += norm[n]; 
00222                 face_count[con[dim*n+d]] += static_cast<DataType>(1.);
00223               }
00224 
00225             for (n=0; n<Npoints; n++) {
00226               point_normal[n] /= face_count[n];
00227               positions[2*n]   = pos[n]-distance*point_normal[n];
00228               positions[2*n+1] = pos[n]+distance*point_normal[n];
00229               centroids[2*n+1] = centroids[2*n] = pos[n];
00230             }
00231 
00232           END_WATCH(FLUID_CPL_ELC_GEOMETRY)
00233           delete [] face_count;
00234           delete [] point_normal;
00235         }
00236         else {
00237           START_WATCH
00238             for (n=0; n<Npoints; n++) {
00239               positions[2*n]   = cent[n]-distance*norm[n];
00240               positions[2*n+1] = cent[n]+distance*norm[n];
00241               centroids[2*n+1] = centroids[2*n] = cent[n];
00242             }
00243           END_WATCH(FLUID_CPL_ELC_GEOMETRY)
00244         }
00245 #ifdef DEBUG_PRINT
00246         ( comm_service::log() << "*** Pressure interpolation()\n" ).flush();
00247 #endif
00248         DataType* press_int = new DataType[2*Npoints];
00249 
00250         START_WATCH
00251           int TimeZero = CurrentTime(base::GH(),0);
00252           int Lev=0, Np=2*Npoints;
00253           bool* lpos = new bool[Np];
00254           bool* lcent = new bool[Np];
00255           _Interpolation->LocalPoints(base::Work(),TimeZero,Lev,Np,positions,lpos);
00256           _Interpolation->LocalPoints(base::Work(),TimeZero,Lev,Np,centroids,lcent);
00257           int nlpos=0, nppos=0; 
00258           for (n=0; n<Np; n++)
00259             if (lcent[n])
00260               if (lpos[n]) nlpos++;
00261               else nppos++;
00262           int lm, pm;
00263           point_type* lpositions = new point_type[nlpos];
00264           point_type* ppositions = new point_type[nppos];
00265           DataType* lpress_int = new DataType[nlpos];
00266           DataType* ppress_int = new DataType[nppos];
00267           for (n=0, lm=0, pm=0; n<Np; n++)
00268             if (lcent[n]) 
00269               if (lpos[n]) lpositions[lm++] = positions[n];
00270               else ppositions[pm++] = positions[n];
00271           _Interpolation->PointsValues(base::Work(),base::t[0],nlpos,lpositions,lpress_int);
00272           _Interpolation->PointsValuesPar(base::Work(),base::t[0],nppos,ppositions,ppress_int);
00273           for (n=0, lm=0, pm=0; n<Np; n++)
00274             if (lcent[n]) 
00275               if (lpos[n]) press_int[n] = lpress_int[lm++];
00276               else press_int[n] = ppress_int[pm++];
00277             else press_int[n] = _Interpolation->ErrorValue();
00278           delete [] lpositions;
00279           delete [] ppositions;
00280           delete [] lpress_int;
00281           delete [] ppress_int;
00282           delete [] lpos;
00283           delete [] lcent;
00284         END_WATCH(FLUID_CPL_INTERPOLATE)
00285           
00286         for (n=0; n<Npoints; n++)
00287           if (press_int[2*n]==_Interpolation->ErrorValue() || 
00288               press_int[2*n+1]==_Interpolation->ErrorValue()) {
00289             press_data[n] = std::numeric_limits<DataType>::max();
00290 #ifdef DEBUG_PRINT
00291             if (press_int[2*n]!=press_int[2*n+1] && base::GFM(CoupleGFM()).Boundary().GetVerbose())
00292               ( comm_service::log() << "*** Thin structure interpolation error for center at " 
00293                 << centroids[2*n] << "   " << positions[2*n] << ": " << press_int[2*n] << "   " 
00294                 << positions[2*n+1] << ": " << press_int[2*n+1] << "\n" ).flush();
00295 #endif
00296           }
00297           else
00298             press_data[n] = press_int[2*n]-press_int[2*n+1];
00299         
00300         delete [] positions;
00301         delete [] centroids;
00302         delete [] press_int;
00303       }
00304       else {
00305         if (SendDataLocation()==1) {
00306           if (distance>1.e-12) {
00307             int Nfaces = _eulComm->getNumberOfFaces();
00308             DataType* face_count = new DataType[Npoints];
00309             point_type* point_normal = new point_type[Npoints];
00310             START_WATCH
00311               for (n=0; n<Npoints; n++) {
00312                 face_count[n] = static_cast<DataType>(0.);
00313                 point_normal[n] = static_cast<point_type>(0.);
00314               }
00315               
00316               for (n=0; n<Nfaces; n++) 
00317                 for (register int d=0; d<dim; d++) {
00318                   point_normal[con[dim*n+d]] += norm[n]; 
00319                   face_count[con[dim*n+d]] += static_cast<DataType>(1.);
00320                 }
00321 
00322               point_type* positions = new point_type[Npoints];
00323               for (n=0; n<Npoints; n++) {
00324                 point_normal[n] /= face_count[n];
00325                 positions[n] = pos[n]+distance*point_normal[n];
00326               }
00327             END_WATCH(FLUID_CPL_ELC_GEOMETRY)
00328             delete [] face_count;
00329             delete [] point_normal;
00330 
00331             START_WATCH
00332               _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,positions,press_data);
00333             END_WATCH(FLUID_CPL_INTERPOLATE)
00334 
00335             delete [] positions;
00336           }
00337           else {
00338             START_WATCH
00339               _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,pos,press_data);
00340             END_WATCH(FLUID_CPL_INTERPOLATE)
00341           }
00342         }
00343         else {
00344           point_type* positions = new point_type[Npoints];
00345           START_WATCH
00346             for (n=0; n<Npoints; n++) { 
00347               positions[n] = cent[n];
00348               if (distance>1.e-12) 
00349                 positions[n] += distance*norm[n];
00350             }   
00351           END_WATCH(FLUID_CPL_ELC_GEOMETRY)
00352         
00353           START_WATCH
00354             _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,positions,press_data);
00355           END_WATCH(FLUID_CPL_INTERPOLATE)
00356 
00357           delete [] positions;
00358         }
00359         for (n=0; n<Npoints; n++) 
00360           if (press_data[n]==_Interpolation->ErrorValue()) 
00361             press_data[n] = std::numeric_limits<DataType>::max();
00362       }
00363       
00364 #ifdef DEBUG_PRINT_ELC
00365       ( comm_service::log() << "*** ELC::SendPressure" ).flush();
00366 #endif
00367     
00368       if (PressureConversionFactor()!=1.) {
00369         for (n=0; n<Npoints; n++) 
00370           if (press_data[n]!=std::numeric_limits<DataType>::max())
00371             press_data[n] *= static_cast<DataType>(PressureConversionFactor());
00372       }
00373 
00374       ModifyELCSendData(_eulComm);
00375 
00376       START_WATCH
00377         _eulComm->sendPressure();
00378         _eulComm->waitForPressure();
00379       END_WATCH(ELC_SENDPRESSURE)
00380 
00381 #ifdef DEBUG_PRINT_ELC
00382       ( comm_service::log() << " done.\n" ).flush();
00383 #endif
00384 #ifdef DEBUG_PRINT_ELC_VALUES
00385       ( comm_service::log() << "*** Values passed to ELC ***\n" ).flush();
00386       for (n=0; n<Npoints; n++) 
00387         ( comm_service::log() << press_data[n] << "\n" ).flush();    
00388 #endif
00389     END_WATCH(FLUID_CPL_SEND_OVERHEAD)
00390   }
00391 
00392   void ExtractBoundaryVelocities(const BBox& bb, const int& nc, 
00393                                  const point_type* xc, point_type* retvel) {
00394     START_WATCH
00395       int nodes = _eulComm->getNumberOfNodes();
00396       const point_type *positions = reinterpret_cast<const point_type *>(_eulComm->getPositionsData());
00397       const point_type *velocities = reinterpret_cast<const point_type *>(_eulComm->getVelocitiesData());
00398 
00399 #ifdef DEBUG_PRINT_ELC_VALUES
00400       ( comm_service::log() << "*** Velocities constructed in " << bb << " ***\n" ).flush();
00401 #endif
00402 
00403       // We are looking here only in a Box, which is not necessarily enough
00404       // temporary solution until CPT is used for that
00405       DCoords lbcorner = base::GH().worldCoords(bb.lower(),bb.stepsize());
00406       DCoords ubcorner = base::GH().worldCoords(bb.upper()+bb.stepsize(),bb.stepsize());
00407       point_type lb, ub;
00408       register int m, n;
00409       for (m=0; m<base::Dim(); m++) {
00410         lb(m) = lbcorner(m); ub(m) = ubcorner(m);
00411       }
00412 
00413       std::vector<int> pos_list;
00414       for (m=0; m<nodes; m++)
00415         if (!(lb>positions[m] || ub<positions[m]))
00416           pos_list.push_back(m);    
00417       int s = pos_list.size();
00418 
00419       if (s>0) {
00420         double min_dist, dist;
00421         register int min_point;
00422         for (n=0; n<nc; n++) {
00423           min_dist = Abs(xc[n]-positions[pos_list[0]]);
00424           min_point = pos_list[0];
00425           for (m=1; m<s; m++) {
00426             dist = Abs(xc[n]-positions[pos_list[m]]);
00427             if (dist<min_dist) {
00428               min_dist = dist;
00429               min_point = pos_list[m];
00430             }
00431           }
00432           retvel[n] = velocities[min_point];      
00433         }
00434       }
00435       else 
00436         for (n=0; n<nc; n++) 
00437           retvel[n] = 0.;
00438    
00439 #ifdef DEBUG_PRINT_ELC_VALUES
00440       for (n=0; n<nc; n++) 
00441         ( comm_service::log() << xc[n] << "   " << retvel[n] << "\n" ).flush();
00442 #endif
00443     END_WATCH(FLUID_CPL_VELOCITY_SEARCH)
00444   }
00445 
00446   virtual void ModifyELCReceiveData(eul_comm_type* eulComm) {}
00447   virtual void ModifyELCSendData(eul_comm_type* eulComm) {}
00448 
00449   inline void SetCoupleGFM(int gfm) { _CoupleGFM = gfm; }
00450   inline const int& CoupleGFM() const { return _CoupleGFM; }
00451   inline void SetSendDataLocation(int bd) { SendDataLocation = bd; }
00452   inline const int& SendDataLocation() const { return _SendDataLocation; }
00453   inline void SetLengthConversionFactor(DataType l) { _LengthConversionFactor = l; }
00454   inline const DataType& LengthConversionFactor() const { return _LengthConversionFactor; }
00455   inline void SetPressureConversionFactor(DataType p) { _PressureConversionFactor = p; }
00456   inline const DataType& PressureConversionFactor() const { return _PressureConversionFactor; }
00457 
00458 protected:
00459   interpolation_type* _Interpolation;
00460   eul_comm_type* _eulComm;
00461   int _CoupleGFM;
00462   int rank;
00463   int _SendDataLocation, _ThinStructure;
00464   DataType _LengthConversionFactor, _PressureConversionFactor;
00465 };
00466 
00467 
00468 template <class VectorType, int dim>
00469 class F77ELCGFMBoundary : public F77GFMBoundary<VectorType,dim> {
00470   typedef typename VectorType::InternalDataType DataType;
00471   typedef F77GFMBoundary<VectorType,dim> base;
00472 public:
00473   typedef typename base::vec_grid_data_type vec_grid_data_type;
00474   typedef typename base::grid_data_type grid_data_type;
00475   typedef typename base::generic_func_type generic_func_type;
00476   typedef typename base::point_type point_type;
00477   typedef AMRELCGFMSolver<VectorType,FixupType,FlagType,dim> amr_elc_solver;
00478 
00479   F77ELCGFMBoundary(generic_func_type bnd, generic_func_type trs, amr_elc_solver& solver) : 
00480     base(bnd,trs), _solver(solver) {
00481     base::SetNAux(base::Dim());
00482   }
00483 
00484   virtual void SetBndryAux(vec_grid_data_type& gdu, grid_data_type& gdphi, const BBox& bb,
00485                            const VectorType* u, DataType* aux, const int& Level, 
00486                            double t, const int& nc, const int* idx, 
00487                            const point_type* xc, const DataType* distance, 
00488                            const point_type* normal) {
00489     _solver.ExtractBoundaryVelocities(gdu.bbox(),nc,xc,reinterpret_cast<point_type *>(aux));
00490   }
00491 
00492 protected:
00493   amr_elc_solver& _solver;
00494 };
00495 
00496 
00497 #endif

Generated on Fri Aug 24 13:00:45 2007 for AMROC Fluid-solver Framework - by  doxygen 1.4.7