Skip to content

Commit d916987

Browse files
asterioukonstantinosasterioukonstantinos
committed
Fully working version for checkpointing when preinlet exists and solidification or viscosity is enabled
1 parent 3189621 commit d916987

File tree

5 files changed

+127
-51
lines changed

5 files changed

+127
-51
lines changed

‎config/constant_defaults.h

+1-1
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
3232
* Which version are we at
3333
*/
3434
#define VERSION_MAJOR 2
35-
#define VERSION_MINOR "6"
35+
#define VERSION_MINOR "7"
3636

3737
// Use for solidifying, performance impact so only enable when using
3838
#ifndef SOLIDIFY_MECHANICS

‎core/hemoCell.cpp

+13-12
Original file line numberDiff line numberDiff line change
@@ -327,7 +327,6 @@ void HemoCell::iterate() {
327327
if(iter %cellfields->particleVelocityUpdateTimescale == 0) {
328328
// #### 3 #### IBM interpolation
329329
cellfields->interpolateFluidVelocity();
330-
331330
// ### 4 ### sync the particles
332331
cellfields->syncEnvelopes();
333332
}
@@ -482,18 +481,20 @@ void HemoCell::initializeLattice(MultiBlockManagement3D const & management) {
482481
totalNodes += preInlet->getNumberOfNodes();
483482
totalNodes += cellsInBoundingBox(management.getBoundingBox());
484483

485-
preInlet->nProcs = global::mpi().getSize()*(preInlet->getNumberOfNodes()/(T)totalNodes);
486-
if (preInlet->nProcs == 0) {
487-
preInlet->nProcs = 1;
484+
try { // Look for block management info in the config file
485+
plint preInlet_pABx = (*cfg)["preInlet"]["parameters"]["pABx"].read<plint>();
486+
plint preInlet_pABy = (*cfg)["preInlet"]["parameters"]["pABy"].read<plint>();
487+
plint preInlet_pABz = (*cfg)["preInlet"]["parameters"]["pABz"].read<plint>();
488+
preInlet->nProcs = preInlet_pABx*preInlet_pABy*preInlet_pABz;
489+
}
490+
catch (const std::invalid_argument& e) {
491+
preInlet->nProcs = global::mpi().getSize()*(preInlet->getNumberOfNodes()/(T)totalNodes);
492+
if (preInlet->nProcs == 0) {
493+
preInlet->nProcs = 1;
494+
}
488495
}
489-
490-
// JON addition: Try to read in number of processors allocated to preinlet from config file.
491-
// TODO: Do we still need this?
492-
// Just continue with value computed above if reading it from XML throws an exception because it does not exist
493-
try { preInlet->nProcs = (*cfg)["preInlet"]["parameters"]["nProcs"].read<int>(); }
494-
catch (const std::invalid_argument& e) {}
495496

496-
int nProcs = global::mpi().getSize()-preInlet->nProcs;
497+
int nProcs = global::mpi().getSize() - preInlet->nProcs;
497498

498499
//Assign processors to PreInlet or Domain
499500
unsigned int currentPreInlet = 0;
@@ -692,4 +693,4 @@ void HemoCell::sanityCheck() {
692693

693694

694695
sanityCheckDone = true;
695-
}
696+
}

‎helper/bindingField.cpp

+18-14
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
2525
#include "hemocell.h"
2626
#include "palabos3D.h"
2727
#include "palabos3D.hh"
28+
#include "preInlet.h"
2829

2930
namespace hemo {
3031
bindingFieldHelper::bindingFieldHelper(HemoCellFields * cellFields_) : cellFields(*cellFields_) {
@@ -36,26 +37,26 @@ namespace hemo {
3637

3738
//Create bindingfield with same properties as fluid field underlying the particleField.
3839
multiBindingField = new plb::MultiScalarField3D<bool>(
39-
MultiBlockManagement3D (
40-
*cellFields.hemocell.lattice->getSparseBlockStructure().clone(),
41-
cellFields.hemocell.lattice->getMultiBlockManagement().getThreadAttribution().clone(),
42-
cellFields.hemocell.lattice->getMultiBlockManagement().getEnvelopeWidth(),
43-
cellFields.hemocell.lattice->getMultiBlockManagement().getRefinementLevel()),
40+
MultiBlockManagement3D (
41+
*cellFields.hemocell.domain_lattice->getSparseBlockStructure().clone(),
42+
cellFields.hemocell.domain_lattice->getMultiBlockManagement().getThreadAttribution().clone(),
43+
cellFields.hemocell.domain_lattice->getMultiBlockManagement().getEnvelopeWidth(),
44+
cellFields.hemocell.domain_lattice->getMultiBlockManagement().getRefinementLevel()),
4445
defaultMultiBlockPolicy3D().getBlockCommunicator(),
4546
defaultMultiBlockPolicy3D().getCombinedStatistics(),
4647
defaultMultiBlockPolicy3D().getMultiScalarAccess<bool>(),
4748
0);
48-
multiBindingField->periodicity().toggle(0,cellFields.hemocell.lattice->periodicity().get(0));
49-
multiBindingField->periodicity().toggle(1,cellFields.hemocell.lattice->periodicity().get(1));
50-
multiBindingField->periodicity().toggle(2,cellFields.hemocell.lattice->periodicity().get(2));
49+
multiBindingField->periodicity().toggle(0,cellFields.hemocell.domain_lattice->periodicity().get(0));
50+
multiBindingField->periodicity().toggle(1,cellFields.hemocell.domain_lattice->periodicity().get(1));
51+
multiBindingField->periodicity().toggle(2,cellFields.hemocell.domain_lattice->periodicity().get(2));
5152

5253
multiBindingField->initialize();
53-
5454
//Make sure each particleField has access to its local scalarField
5555
for (const plint & bId : multiBindingField->getLocalInfo().getBlocks()) {
56-
HemoCellParticleField & pf = cellFields.immersedParticles->getComponent(bId);
56+
HemoCellParticleField & pf = cellFields.domain_immersedParticles->getComponent(bId);
5757
pf.bindingField = &multiBindingField->getComponent(bId);
5858
}
59+
5960
}
6061

6162
bindingFieldHelper::~bindingFieldHelper() {
@@ -67,15 +68,17 @@ namespace hemo {
6768
pcout << "(BindingField) Checkpoint called while global.enableSolidifyMechanics is not enabled, still checkpointing but this should not happen" << endl;
6869
return;
6970
}
71+
7072
std::string & outDir = hemo::global.checkpointDirectory;
7173
mkpath(outDir.c_str(), 0777);
72-
74+
7375
if (global::mpi().isMainProcessor()) {
7476
renameFileToDotOld(outDir + "bindingSites.dat");
7577
renameFileToDotOld(outDir + "bindingSites.plb");
7678
}
77-
79+
7880
plb::parallelIO::save(*multiBindingField, outDir + "bindingSites", true);
81+
7982
}
8083

8184
void bindingFieldHelper::restore(HemoCellFields & cellFields) {
@@ -119,8 +122,8 @@ namespace hemo {
119122
}
120123

121124
void bindingFieldHelper::refillBindingSites() {
122-
for (const plint & bId : cellFields.immersedParticles->getLocalInfo().getBlocks()) {
123-
HemoCellParticleField & pf = cellFields.immersedParticles->getComponent(bId);
125+
for (const plint & bId : cellFields.domain_immersedParticles->getLocalInfo().getBlocks()) {
126+
HemoCellParticleField & pf = cellFields.domain_immersedParticles->getComponent(bId);
124127
ScalarField3D<bool> & bf = *pf.bindingField;
125128
Box3D domain = bf.getBoundingBox();
126129
for (int x = domain.x0; x <= domain.x1 ; x++) {
@@ -134,4 +137,5 @@ namespace hemo {
134137
}
135138
}
136139
}
140+
137141
}

‎helper/interiorViscosity.cpp

+89-21
Original file line numberDiff line numberDiff line change
@@ -28,28 +28,59 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
2828

2929
namespace hemo {
3030
InteriorViscosityHelper::InteriorViscosityHelper(HemoCellFields & cellFields_) : cellFields(cellFields_) {
31-
//Create bindingfield with same properties as fluid field underlying the particleField.
32-
multiInteriorViscosityField = new plb::MultiScalarField3D<T>(
31+
//Create viscosity field with same properties as fluid field underlying the particleField.
32+
if(cellFields.hemocell.preInlet){
33+
preinlet_multiInteriorViscosityField = new plb::MultiScalarField3D<T>(
3334
MultiBlockManagement3D (
34-
*cellFields.hemocell.lattice->getSparseBlockStructure().clone(),
35-
cellFields.hemocell.lattice->getMultiBlockManagement().getThreadAttribution().clone(),
36-
cellFields.hemocell.lattice->getMultiBlockManagement().getEnvelopeWidth(),
37-
cellFields.hemocell.lattice->getMultiBlockManagement().getRefinementLevel()),
35+
*cellFields.hemocell.preinlet_lattice->getSparseBlockStructure().clone(),
36+
cellFields.hemocell.preinlet_lattice->getMultiBlockManagement().getThreadAttribution().clone(),
37+
cellFields.hemocell.preinlet_lattice->getMultiBlockManagement().getEnvelopeWidth(),
38+
cellFields.hemocell.preinlet_lattice->getMultiBlockManagement().getRefinementLevel()),
3839
defaultMultiBlockPolicy3D().getBlockCommunicator(),
3940
defaultMultiBlockPolicy3D().getCombinedStatistics(),
4041
defaultMultiBlockPolicy3D().getMultiScalarAccess<T>(),
4142
0);
42-
multiInteriorViscosityField->periodicity().toggle(0,cellFields.hemocell.lattice->periodicity().get(0));
43-
multiInteriorViscosityField->periodicity().toggle(1,cellFields.hemocell.lattice->periodicity().get(1));
44-
multiInteriorViscosityField->periodicity().toggle(2,cellFields.hemocell.lattice->periodicity().get(2));
43+
preinlet_multiInteriorViscosityField->periodicity().toggle(0,cellFields.hemocell.preinlet_lattice->periodicity().get(0));
44+
preinlet_multiInteriorViscosityField->periodicity().toggle(1,cellFields.hemocell.preinlet_lattice->periodicity().get(1));
45+
preinlet_multiInteriorViscosityField->periodicity().toggle(2,cellFields.hemocell.preinlet_lattice->periodicity().get(2));
4546

46-
multiInteriorViscosityField->initialize();
47+
preinlet_multiInteriorViscosityField->initialize();
48+
49+
//Make sure each particleField has access to its local scalarField
50+
for (const plint & bId : preinlet_multiInteriorViscosityField->getLocalInfo().getBlocks()) {
51+
HemoCellParticleField & pf = cellFields.preinlet_immersedParticles->getComponent(bId);
52+
pf.interiorViscosityField = &preinlet_multiInteriorViscosityField->getComponent(bId);
53+
}
54+
}
55+
domain_multiInteriorViscosityField = new plb::MultiScalarField3D<T>(
56+
MultiBlockManagement3D (
57+
*cellFields.hemocell.domain_lattice->getSparseBlockStructure().clone(),
58+
cellFields.hemocell.domain_lattice->getMultiBlockManagement().getThreadAttribution().clone(),
59+
cellFields.hemocell.domain_lattice->getMultiBlockManagement().getEnvelopeWidth(),
60+
cellFields.hemocell.domain_lattice->getMultiBlockManagement().getRefinementLevel()),
61+
defaultMultiBlockPolicy3D().getBlockCommunicator(),
62+
defaultMultiBlockPolicy3D().getCombinedStatistics(),
63+
defaultMultiBlockPolicy3D().getMultiScalarAccess<T>(),
64+
0);
65+
domain_multiInteriorViscosityField->periodicity().toggle(0,cellFields.hemocell.domain_lattice->periodicity().get(0));
66+
domain_multiInteriorViscosityField->periodicity().toggle(1,cellFields.hemocell.domain_lattice->periodicity().get(1));
67+
domain_multiInteriorViscosityField->periodicity().toggle(2,cellFields.hemocell.domain_lattice->periodicity().get(2));
68+
69+
domain_multiInteriorViscosityField->initialize();
4770

4871
//Make sure each particleField has access to its local scalarField
49-
for (const plint & bId : multiInteriorViscosityField->getLocalInfo().getBlocks()) {
50-
HemoCellParticleField & pf = cellFields.immersedParticles->getComponent(bId);
51-
pf.interiorViscosityField = &multiInteriorViscosityField->getComponent(bId);
72+
for (const plint & bId : domain_multiInteriorViscosityField->getLocalInfo().getBlocks()) {
73+
HemoCellParticleField & pf = cellFields.domain_immersedParticles->getComponent(bId);
74+
pf.interiorViscosityField = &domain_multiInteriorViscosityField->getComponent(bId);
75+
}
76+
77+
if(cellFields.hemocell.partOfpreInlet){
78+
multiInteriorViscosityField = preinlet_multiInteriorViscosityField;
79+
}
80+
else{
81+
multiInteriorViscosityField = domain_multiInteriorViscosityField;
5282
}
83+
5384
}
5485

5586
InteriorViscosityHelper::~InteriorViscosityHelper() {
@@ -67,9 +98,15 @@ namespace hemo {
6798
if (global::mpi().isMainProcessor()) {
6899
renameFileToDotOld(outDir + "internalViscosity.dat");
69100
renameFileToDotOld(outDir + "internalViscosity.plb");
101+
if(cellFields.hemocell.preInlet){
102+
renameFileToDotOld(outDir + "PRE_internalViscosity.dat");
103+
renameFileToDotOld(outDir + "PRE_internalViscosity.plb");
104+
}
70105
}
71-
72-
plb::parallelIO::save(*multiInteriorViscosityField, outDir + "internalViscosity", true);
106+
if(cellFields.hemocell.preInlet){
107+
plb::parallelIO::save(*preinlet_multiInteriorViscosityField, outDir + "PRE_internalViscosity", true);
108+
}
109+
plb::parallelIO::save(*domain_multiInteriorViscosityField, outDir + "internalViscosity", true);
73110
}
74111

75112
void InteriorViscosityHelper::restore(HemoCellFields & cellFields) {
@@ -84,8 +121,18 @@ namespace hemo {
84121
pcout << "(internalViscosityField) Error restoring internalViscosity fields from checkpoint, they do not seem to exist" << endl;
85122
exit(1);
86123
}
87-
88-
plb::parallelIO::load(outDir + "internalViscosity",*get(cellFields).multiInteriorViscosityField,true);
124+
if(cellFields.hemocell.preInlet){
125+
std::string file_dat = outDir + "PRE_internalViscosity.dat";
126+
std::string file_plb = outDir + "PRE_internalViscosity.plb";
127+
if(!(file_exists(file_dat) && file_exists(file_plb))) {
128+
pcout << "(PRE_internalViscosityField) Error restoring PRE_internalViscosity fields from checkpoint, they do not seem to exist" << endl;
129+
exit(1);
130+
}
131+
}
132+
plb::parallelIO::load(outDir + "internalViscosity",*get(cellFields).domain_multiInteriorViscosityField,true);
133+
if(cellFields.hemocell.preInlet){
134+
plb::parallelIO::load(outDir + "PRE_internalViscosity",*get(cellFields).preinlet_multiInteriorViscosityField,true);
135+
}
89136
get(cellFields).refillBindingSites();
90137
}
91138

@@ -119,8 +166,29 @@ namespace hemo {
119166
}
120167

121168
void InteriorViscosityHelper::refillBindingSites() {
122-
for (const plint & bId : cellFields.immersedParticles->getLocalInfo().getBlocks()) {
123-
HemoCellParticleField & pf = cellFields.immersedParticles->getComponent(bId);
169+
if(cellFields.hemocell.preInlet){
170+
for (const plint & bId : cellFields.preinlet_immersedParticles->getLocalInfo().getBlocks()) {
171+
HemoCellParticleField & pf = cellFields.preinlet_immersedParticles->getComponent(bId);
172+
ScalarField3D<T> & bf = *pf.interiorViscosityField;
173+
Box3D domain = bf.getBoundingBox();
174+
for (int x = domain.x0; x <= domain.x1 ; x++) {
175+
for (int y = domain.y0; y <= domain.y1; y++) {
176+
for (int z = domain.z0; z <= domain.z1; z++) {
177+
if(bf.get(x,y,z)) {
178+
pf.internalPoints.insert({x,y,z});
179+
180+
//WARNING this _can_ memory leak, so you should be ok if this is only called one time (from checkpointing)
181+
plb::Dynamics<T,DESCRIPTOR>* dynamic = cellFields.hemocell.preinlet_lattice->getBackgroundDynamics().clone();
182+
dynamic->setOmega(1.0/bf.get(x,y,z));
183+
pf.atomicLattice->get(x,y,z).attributeDynamics(dynamic);
184+
}
185+
}
186+
}
187+
}
188+
}
189+
}
190+
for (const plint & bId : cellFields.domain_immersedParticles->getLocalInfo().getBlocks()) {
191+
HemoCellParticleField & pf = cellFields.domain_immersedParticles->getComponent(bId);
124192
ScalarField3D<T> & bf = *pf.interiorViscosityField;
125193
Box3D domain = bf.getBoundingBox();
126194
for (int x = domain.x0; x <= domain.x1 ; x++) {
@@ -130,7 +198,7 @@ namespace hemo {
130198
pf.internalPoints.insert({x,y,z});
131199

132200
//WARNING this _can_ memory leak, so you should be ok if this is only called one time (from checkpointing)
133-
plb::Dynamics<T,DESCRIPTOR>* dynamic = cellFields.lattice->getBackgroundDynamics().clone();
201+
plb::Dynamics<T,DESCRIPTOR>* dynamic = cellFields.hemocell.domain_lattice->getBackgroundDynamics().clone();
134202
dynamic->setOmega(1.0/bf.get(x,y,z));
135203
pf.atomicLattice->get(x,y,z).attributeDynamics(dynamic);
136204
}
@@ -139,4 +207,4 @@ namespace hemo {
139207
}
140208
}
141209
}
142-
}
210+
}

‎helper/interiorViscosity.h

+6-3
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,8 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
3030
namespace hemo {
3131
class InteriorViscosityHelper {
3232
public:
33-
static InteriorViscosityHelper& get(HemoCellFields & cellFields) {
33+
static InteriorViscosityHelper& get(HemoCellFields & cellFields) {
34+
3435
static InteriorViscosityHelper instance(cellFields);
3536
return instance;
3637
}
@@ -47,7 +48,9 @@ namespace hemo {
4748

4849
private:
4950
HemoCellFields & cellFields;
50-
plb::MultiScalarField3D<T> * multiInteriorViscosityField = 0;
51+
52+
plb::MultiScalarField3D<T> *multiInteriorViscosityField = nullptr,
53+
*preinlet_multiInteriorViscosityField = nullptr, *domain_multiInteriorViscosityField = nullptr;
5154

5255
InteriorViscosityHelper(HemoCellFields & cellFields);
5356
~InteriorViscosityHelper();
@@ -60,4 +63,4 @@ namespace hemo {
6063
void operator=(InteriorViscosityHelper const&) = delete;
6164
};
6265
}
63-
#endif
66+
#endif

0 commit comments

Comments
 (0)