00001
00002
00003
00004
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
00404
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