/*?This?file?is?part?of?the?Palabos?library.
?*
?*?Copyright?（C）?2011-2017?FlowKit?Sarl
?*?Route?d‘Oron?2
?*?1010?Lausanne?Switzerland
?*?E-mail?contact:?contact@flowkit.com
?*
?*?The?most?recent?release?of?Palabos?can?be?downloaded?at?
?*?
?*
?*?The?library?Palabos?is?free?software:?you?can?redistribute?it?and/or
?*?modify?it?under?the?terms?of?the?GNU?Affero?General?Public?License?as
?*?published?by?the?Free?Software?Foundation?either?version?3?of?the
?*?License?or?（at?your?option）?any?later?version.
?*
?*?The?library?is?distributed?in?the?hope?that?it?will?be?useful
?*?but?WITHOUT?ANY?WARRANTY;?without?even?the?implied?warranty?of
?*?MERCHANTABILITY?or?FITNESS?FOR?A?PARTICULAR?PURPOSE.??See?the
?*?GNU?Affero?General?Public?License?for?more?details.
?*
?*?You?should?have?received?a?copy?of?the?GNU?Affero?General?Public?License
?*?along?with?this?program.??If?not?see?.
*/

/*?The?breaking?dam?free?surface?problem.?This?code?demonstrates?the?basic?usage?of?the
?*?free?surface?module?in?Palabos.?Surface?tension?and?contact?angles?are?optional.?
?*/

#include?“palabos3D.h“
#include?“palabos3D.hh“

using?namespace?plb;

#define?DEscriptOR?descriptors::ForcedD3Q19Descriptor

typedef?double?T;


//?Smagorinsky?constant?for?LES?model.
const?T?cSmago?=?0.14;

//?Physical?dimensions?of?the?system?（in?meters）.
const?T?lx?=?3.22;
const?T?ly?=?1.0;
const?T?lz?=?1.0;

const?T?rhoEmpty?=?T（1）;
????
plint?writeImagesIter???=?100;
plint?getStatisticsIter?=?20;

plint?maxIter;
plint?N;
plint?nx?ny?nz;
T?delta_t?delta_x;
Array?externalForce;
T?nuPhys?nuLB?tau?omega?Bo?surfaceTensionLB?contactAngle;

std::string?outDir;
plint?obstacleCenterXYplane?obstacleLength?obstacleWidth?obstacleHeight?beginWaterReservoir?waterReservoirHeight;
plint?waterLevelOne?waterLevelTwo?waterLevelThree?waterLevelFour;

void?setupParameters（）?{
????delta_x?=?lz?/?N;
????nx?=?util::roundToInt（lx?/?delta_x）;
????ny?=?util::roundToInt（ly?/?delta_x）;
????nz?=?util::roundToInt（lz?/?delta_x）;

????//?Gravity?in?lattice?units.
????T?gLB?=?9.8?*?delta_t?*?delta_t/delta_x;
????externalForce?=?Array（0.?0.?-gLB）;
????tau????????????=?（nuPhys*DEscriptOR::invCs2*delta_t）/（delta_x*delta_x）?+?0.5;
????omega??????????=?1./tau;????
????nuLB???????????=?（tau-0.5）*DEscriptOR::cs2;?//?Viscosity?in?lattice?units.
????
????surfaceTensionLB?=?rhoEmpty?*?gLB?*?N?*?N?/?Bo;

????obstacleCenterXYplane?=?util::roundToInt（0.744*N）;
????obstacleLength????????=?util::roundToInt（0.403*N）;
????obstacleWidth?????????=?util::roundToInt（0.161*N）;
????obstacleHeight????????=?util::roundToInt（0.161*N）;
????beginWaterReservoir???=?util::roundToInt（（0.744+1.248）*N）;
????waterReservoirHeight??=?util::roundToInt（0.55*N）;
????
????waterLevelOne???=?util::roundToInt（0.496*N）;
????waterLevelTwo???=?util::roundToInt（2.*0.496*N）;
????waterLevelThree?=?util::roundToInt（3.*0.496*N）;
????waterLevelFour??=?util::roundToInt（（3.*0.496?+?1.150）*N