/** * @file tempgrid.cpp * @author Julian Iseringhausen * @date Mon Apr 18 12:55:05 2011 * * @brief VMG::TempGrid * */ #ifdef HAVE_CONFIG_H #include #endif #include "base/discretization.hpp" #include "base/stencil.hpp" #include "comm/comm.hpp" #include "grid/grid_index_translations.hpp" #include "grid/tempgrid.hpp" #include "interface/interface.hpp" #include "mg.hpp" using namespace VMG; void TempGrid::SetProperties(const Grid& rhs) { local = rhs.Local(); global = rhs.Global(); extent = rhs.Extent(); iterators.SetSubgrids(rhs.Local()); level = rhs.Level(); Allocate(); } void TempGrid::SetProperties(const GlobalIndices& global_, const LocalIndices& local_, const SpatialExtent& extent_) { local = local_; global = global_; extent = extent_; iterators.SetSubgrids(local_); Allocate(); } void TempGrid::SetPropertiesToGlobalCoarseGrid() { Interface* interface = MG::GetInterface(); global = interface->Global().back(); extent = interface->Extent().back(); for (int i=0; i<3; ++i) { if (interface->BoundaryConditions()[i] == Dirichlet || interface->BoundaryConditions()[i] == Open) { local.Size()[i] = global.SizeGlobal()[i] - 2; local.SizeTotal()[i] = global.SizeGlobal()[i]; local.Begin()[i] = 1; local.End()[i] = local.SizeTotal()[i] - 1; local.HaloBegin1()[i] = 0; local.HaloEnd1()[i] = 0; local.HaloBegin2()[i] = 0; local.HaloEnd2()[i] = 0; local.BoundaryBegin1()[i] = 0; local.BoundaryEnd1()[i] = 1; local.BoundaryBegin2()[i] = local.SizeTotal()[i] - 1; local.BoundaryEnd2()[i] = local.SizeTotal()[i]; }else if (interface->BoundaryConditions()[i] == Periodic) { local.Size()[i] = global.SizeGlobal()[i]; local.SizeTotal()[i] = global.SizeGlobal()[i]; local.Begin()[i] = 0; local.End()[i] = global.SizeGlobal()[i]; local.HaloBegin1()[i] = 0; local.HaloEnd1()[i] = 0; local.HaloBegin2()[i] = 0; local.HaloEnd2()[i] = 0; local.BoundaryBegin1()[i] = 0; local.BoundaryEnd1()[i] = 0; local.BoundaryBegin2()[i] = 0; local.BoundaryEnd2()[i] = 0; } } level = interface->MinLevel(); iterators.SetSubgrids(local); Allocate(); } void TempGrid::SetPropertiesToFiner(const Grid& grid, const Boundary& boundary) { assert(grid.Father() != NULL); assert(grid.Level() < grid.Father()->MaxLevel()); const Index off = GridIndexTranslations::EndOffset(boundary); this->Global().BeginFinest() = grid.Global().BeginFinest(); this->Global().EndFinest() = grid.Global().EndFinest(); this->Global().SizeFinest() = grid.Global().SizeFinest(); this->Global().SizeGlobal() = 2 * (grid.Global().SizeGlobal()-off) + off; this->Global().BeginLocal() = grid.Global().BeginLocal(); this->Global().EndLocal() = grid.Global().EndLocal(); GridIndexTranslations::CoarseToFine(this->Global().BeginLocal(), this->Global().EndLocal(), this->Global().SizeGlobal()); this->Global().SizeLocal() = this->Global().EndLocal() - this->Global().BeginLocal(); this->Global().BoundaryType() = (*(grid.Father()))(grid.Level()+1).Global().BoundaryType(); this->Local().Size() = this->Global().SizeLocal(); this->Local().SizeTotal() = this->Global().SizeLocal(); this->Local().Begin() = 0; this->Local().End() = this->Global().SizeLocal(); for (int i=0; i<3; ++i) { if (grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i] > 0) { this->Local().SizeTotal()[i] += 1; this->Local().Begin()[i] += 1; this->Local().End()[i] += 1; this->Local().HaloBegin1()[i] = 0; this->Local().HaloEnd1()[i] = 1; }else { this->Local().HaloBegin1()[i] = 0; this->Local().HaloEnd1()[i] = 0; } if (grid.Local().BoundaryEnd1()[i] - grid.Local().BoundaryBegin1()[i] > 0) { this->Local().Size()[i] -= 1; this->Local().Begin()[i] += 1; this->Local().BoundaryBegin1()[i] = 0; this->Local().BoundaryEnd1()[i] = 1; }else { this->Local().BoundaryBegin1()[i] = 0; this->Local().BoundaryEnd1()[i] = 0; } if (grid.Local().HaloEnd2()[i] - grid.Local().HaloBegin2()[i] > 0) { this->Local().SizeTotal()[i] += 1; this->Local().HaloBegin2()[i] = this->Local().End()[i]; this->Local().HaloEnd2()[i] = this->Local().End()[i] + 1; }else { this->Local().HaloBegin2()[i] = 0; this->Local().HaloEnd2()[i] = 0; } if (grid.Local().BoundaryEnd2()[i] - grid.Local().BoundaryBegin2()[i] > 0) { this->Local().Size()[i] -= 1; this->Local().End()[i] -= 1; this->Local().BoundaryBegin2()[i] = this->Local().End()[i]; this->Local().BoundaryEnd2()[i] = this->Local().End()[i]+1; }else { this->Local().BoundaryBegin2()[i] = 0; this->Local().BoundaryEnd2()[i] = 0; } } this->Local().AlignmentBegin() = this->Local().Begin(); this->Local().AlignmentEnd() = this->local.End(); this->Extent().Size() = grid.Extent().Size(); this->Extent().SizeFactor() = grid.Extent().SizeFactor(); this->Extent().Begin() = grid.Extent().Begin(); this->Extent().End() = grid.Extent().End(); this->Extent().MeshWidth() = 0.5 * grid.Extent().MeshWidth(); this->level = grid.Level() + 1; iterators.SetSubgrids(local); Allocate(); } void TempGrid::SetPropertiesToCoarser(const Grid& grid, const Boundary& boundary) { assert(grid.Father() != NULL); assert(grid.Level() > grid.Father()->MinLevel()); const Index off = GridIndexTranslations::EndOffset(boundary); this->level = grid.Level() - 1; this->Global().BeginFinest() = grid.Global().BeginFinest(); this->Global().EndFinest() = grid.Global().EndFinest(); this->Global().SizeFinest() = grid.Global().SizeFinest(); this->Global().BeginLocal() = grid.Global().BeginLocal(); this->Global().EndLocal() = grid.Global().EndLocal(); GridIndexTranslations::FineToCoarse(this->Global().BeginLocal(), this->Global().EndLocal()); this->Global().SizeLocal() = this->Global().EndLocal() - this->Global().BeginLocal(); this->Global().SizeGlobal() = (grid.Global().SizeGlobal()-off)/2 + off; this->Global().BoundaryType() = (*(grid.Father()))(grid.Level()-1).Global().BoundaryType(); this->Local().SizeTotal() = this->Global().SizeLocal(); this->Local().Size() = this->Global().SizeLocal(); this->Local().Begin() = 0; this->Local().End() = this->Global().SizeLocal(); for (int i=0; i<3; ++i) { if (grid.Local().HaloEnd1()[i] - grid.Local().HaloBegin1()[i] > 0) { this->Local().SizeTotal()[i] += 1; this->Local().Begin()[i] += 1; this->Local().End()[i] += 1; this->Local().HaloBegin1()[i] = 0; this->Local().HaloEnd1()[i] = 1; }else { this->Local().HaloBegin1()[i] = 0; this->Local().HaloEnd1()[i] = 0; } if (grid.Local().BoundaryEnd1()[i] - grid.Local().BoundaryBegin1()[i] > 0) { this->Local().Size()[i] -= 1; this->Local().Begin()[i] += 1; this->Local().BoundaryBegin1()[i] = 0; this->Local().BoundaryEnd1()[i] = 1; }else { this->Local().BoundaryBegin1()[i] = 0; this->Local().BoundaryEnd1()[i] = 0; } if (grid.Local().HaloEnd2()[i] - grid.Local().HaloBegin2()[i] > 0) { this->Local().SizeTotal()[i] += 1; this->Local().HaloBegin2()[i] = this->Local().End()[i]; this->Local().HaloEnd2()[i] = this->Local().End()[i] + 1; }else { this->Local().HaloBegin2()[i] = 0; this->Local().HaloEnd2()[i] = 0; } if (grid.Local().BoundaryEnd2()[i] - grid.Local().BoundaryBegin2()[i] > 0) { this->Local().Size()[i] -= 1; this->Local().End()[i] -= 1; this->Local().BoundaryBegin2()[i] = this->Local().End()[i]; this->Local().BoundaryEnd2()[i] = this->Local().End()[i]+1; }else { this->Local().BoundaryBegin2()[i] = 0; this->Local().BoundaryEnd2()[i] = 0; } } this->Local().AlignmentBegin() = this->Local().Begin(); this->Local().AlignmentEnd() = this->local.End(); this->Extent().Size() = grid.Extent().Size(); this->Extent().SizeFactor() = grid.Extent().SizeFactor(); this->Extent().Begin() = grid.Extent().Begin(); this->Extent().End() = grid.Extent().End(); this->Extent().MeshWidth() = 2.0 * grid.Extent().MeshWidth(); iterators.SetSubgrids(local); Allocate(); } void TempGrid::ImportFromResidual(Grid& sol, Grid& rhs) { #ifdef DEBUG this->IsCompatible(sol); this->IsCompatible(rhs); #endif Grid::iterator grid_iter; Stencil::iterator stencil_iter; vmg_float residual; Index i; const vmg_float prefactor = MG::GetDiscretization()->OperatorPrefactor(sol); const Stencil& A = MG::GetDiscretization()->GetStencil(); this->Clear(); MG::GetComm()->CommToGhosts(sol); for (grid_iter=Iterators().Local().Begin(); grid_iter!=Iterators().Local().End(); ++grid_iter) { residual = rhs.GetVal(*grid_iter) - prefactor * A.GetDiag() * sol.GetVal(*grid_iter); for (stencil_iter=A.begin(); stencil_iter!=A.end(); ++stencil_iter) residual -= prefactor * stencil_iter->Val() * sol.GetVal(*grid_iter + stencil_iter->Disp()); (*this)(*grid_iter) = residual; } MG::GetComm()->CommToGhosts(*this); } void TempGrid::Allocate() { const int size = local.SizeTotal().Product(); if (size > size_max) { size_max = size; delete [] grid; grid = new vmg_float[size]; } } TempGrid::TempGrid() { size_max = 0; } TempGrid::~TempGrid() { }