Skip to content

feat: add fixed_axis keyword to enable 2d pressure coupling #526

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ape33 wants to merge 1 commit into
base: dev
Choose a base branch
from
Open

feat: add fixed_axis keyword to enable 2d pressure coupling #526

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions changes/user/enhancement.2d-pressure-couling.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
- Enable isotropic, anisotropic and full-anisotropic 2d pressure coupling by adding the "fixed_axis" keyword
1 change: 1 addition & 0 deletions include/input/inputFileParser/manostatInputParser.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@ namespace input
);

void parseIsotropy(const std::vector<std::string> &, const size_t);
void parseFixedAxis(const std::vector<std::string> &, const size_t);
};

} // namespace input
Expand Down
17 changes: 12 additions & 5 deletions include/manostat/berendsenManostat.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,18 @@ namespace manostat
class BerendsenManostat : public Manostat
{
protected:
double _tau;
double _compressibility;
double _dt;
double _tau;
double _compressibility;
double _dt;
settings::FixedAxis _fixedAxis;

public:
explicit BerendsenManostat(const double, const double, const double);
explicit BerendsenManostat(
const double,
const double,
const double,
const settings::FixedAxis
);

void applyManostat(
simulationBox::SimulationBox &simBox,
Expand Down Expand Up @@ -83,7 +89,8 @@ namespace manostat
const double,
const double,
const size_t,
const std::vector<size_t> &
const std::vector<size_t> &,
const settings::FixedAxis
);

[[nodiscard]] linearAlgebra::tensor3D calculateMu() const override;
Expand Down
17 changes: 10 additions & 7 deletions include/manostat/stochasticRescalingManostat.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,16 +40,18 @@ namespace manostat
protected:
randomNumberGenerator::RandomNumberGenerator _randomNumberGenerator{};

double _tau;
double _compressibility;
double _dt;
double _tau;
double _compressibility;
double _dt;
settings::FixedAxis _fixedAxis;

public:
StochasticRescalingManostat() = default;
explicit StochasticRescalingManostat(
const double targetPressure,
const double tau,
const double compressibility
const double targetPressure,
const double tau,
const double compressibility,
const settings::FixedAxis fixedAxis
);
~StochasticRescalingManostat() override = default;

Expand Down Expand Up @@ -97,7 +99,8 @@ namespace manostat
const double tau,
const double compressibility,
const size_t anisotropicAxis,
const std::vector<size_t> &isotropicAxes
const std::vector<size_t> &isotropicAxes,
const settings::FixedAxis fixedAxis
);

[[nodiscard]]
Expand Down
13 changes: 13 additions & 0 deletions include/settings/manostatSettings.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,14 @@ namespace settings
FULL_ANISOTROPIC
};

enum class FixedAxis
{
NONE,
X,
Y,
Z
};

[[nodiscard]] std::string string(const ManostatType &manostatType);
[[nodiscard]] std::string string(const Isotropy &isotropy);

Expand All @@ -74,6 +82,7 @@ namespace settings
private:
static inline ManostatType _manostatType = ManostatType::NONE;
static inline Isotropy _isotropy = Isotropy::ISOTROPIC;
static inline FixedAxis _fixedAxis = FixedAxis::NONE;

static inline double _targetPressure;

Expand All @@ -99,6 +108,9 @@ namespace settings
static void setIsotropy(const std::string_view &isotropy);
static void setIsotropy(const Isotropy &isotropy);

static void setFixedAxis(const std::string_view &fixedAxis);
static void setFixedAxis(const FixedAxis &fixedAxis);

static void setTargetPressure(const double targetPressure);
static void setTauManostat(const double tauManostat);
static void setCompressibility(const double compressibility);
Expand All @@ -113,6 +125,7 @@ namespace settings

[[nodiscard]] static ManostatType getManostatType();
[[nodiscard]] static Isotropy getIsotropy();
[[nodiscard]] static FixedAxis getFixedAxis();
[[nodiscard]] static double getTargetPressure();
[[nodiscard]] static double getTauManostat();
[[nodiscard]] static double getCompressibility();
Expand Down
2 changes: 2 additions & 0 deletions include/setup/manostatSetup.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,8 @@ namespace setup
void setupBerendsenManostat();
void setupStochasticRescalingManostat();

void validateIsotropyFixedAxisCombination() const;

void writeSetupInfo() const;
void writeManostatSelection() const;
void writeBerendsenSetup() const;
Expand Down
42 changes: 42 additions & 0 deletions src/input/inputFileParser/manostatInputParser.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,12 @@ ManostatInputParser::ManostatInputParser(Engine &engine)
bindMember(&ManostatInputParser::parseIsotropy, this),
false
);

addKeyword(
std::string("fixed_axis"),
bindMember(&ManostatInputParser::parseFixedAxis, this),
false
);
}

/**
Expand Down Expand Up @@ -288,3 +294,39 @@ void ManostatInputParser::parseIsotropy(
);
}
}

void ManostatInputParser::parseFixedAxis(
const std::vector<std::string> &lineElements,
const size_t lineNumber
)
{
checkCommand(lineElements, lineNumber);

const auto fixed_axis = toLowerAndReplaceDashesCopy(lineElements[2]);

using enum FixedAxis;

if (fixed_axis == "none")
ManostatSettings::setFixedAxis(NONE);

else if (fixed_axis == "x")
ManostatSettings::setFixedAxis(X);

else if (fixed_axis == "y")
ManostatSettings::setFixedAxis(Y);

else if (fixed_axis == "z")
ManostatSettings::setFixedAxis(Z);

else
{
throw InputFileException(
std::format(
"Invalid fixed_axis \"{}\" at line {} in input file.\n"
"Possible options are: none, x, y, z",
lineElements[2],
lineNumber
)
);
}
}
92 changes: 83 additions & 9 deletions src/manostat/berendsenManostat.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,16 +45,19 @@ using namespace physicalData;
* @param targetPressure
* @param tau
* @param compressibility
* @param fixedAxis
*/
BerendsenManostat::BerendsenManostat(
const double targetPressure,
const double tau,
const double compressibility
const double targetPressure,
const double tau,
const double compressibility,
const FixedAxis fixedAxis
)
: Manostat(targetPressure),
_tau(tau),
_compressibility(compressibility),
_dt(TimingsSettings::getTimeStep())
_dt(TimingsSettings::getTimeStep()),
_fixedAxis(fixedAxis)
{
}

Expand All @@ -64,15 +67,19 @@ BerendsenManostat::BerendsenManostat(
* @param targetPressure
* @param tau
* @param compressibility
* @param anisotropicAxis
* @param isotropicAxes
* @param fixedAxis
*/
SemiIsotropicBerendsenManostat::SemiIsotropicBerendsenManostat(
const double targetPressure,
const double tau,
const double compressibility,
const size_t anisotropicAxis,
const std::vector<size_t> &isotropicAxes
const std::vector<size_t> &isotropicAxes,
const FixedAxis fixedAxis
)
: BerendsenManostat(targetPressure, tau, compressibility),
: BerendsenManostat(targetPressure, tau, compressibility, fixedAxis),
_2DAnisotropicAxis(anisotropicAxis),
_2DIsotropicAxes(isotropicAxes)
{
Expand Down Expand Up @@ -118,14 +125,47 @@ void BerendsenManostat::applyManostat(
/**
* @brief calculate mu as scaling factor for Berendsen manostat (isotropic)
*
* @details If a fixed axis is specified, that axis is not scaled (mu = 1.0)
* and the remaining axes are scaled isotropically
*
* @return tensor3D
*/
tensor3D BerendsenManostat::calculateMu() const
{
const auto p = trace(_pressureTensor) / 3.0;
using enum FixedAxis;

const auto preFactor = _compressibility * _dt / _tau;

return diagonalMatrix(::cbrt(1.0 - preFactor * (_targetPressure - p)));
// 2D pressure coupling
if (_fixedAxis != NONE)
{
const auto fixedAxisIndex = static_cast<size_t>(_fixedAxis) - 1;
const auto p_xyz = diagonal(_pressureTensor);

// Calculate average pressure of non-fixed axes
double p_avg = 0.0;
for (size_t i = 0; i < 3; ++i)
if (i != fixedAxisIndex)
p_avg += p_xyz[i];

p_avg /= 2.0;

// Scale factor for non-fixed axes
const auto mu_2D = ::sqrt(1.0 - preFactor * (_targetPressure - p_avg));

Vec3D mu = {1.0, 1.0, 1.0};
for (size_t i = 0; i < 3; ++i)
if (i != fixedAxisIndex)
mu[i] = mu_2D;

return diagonalMatrix(mu);
}

// 3D pressure coupling
const auto p = trace(_pressureTensor) / 3.0;
const auto mu_scalar = ::cbrt(1.0 - preFactor * (_targetPressure - p));

return diagonalMatrix(mu_scalar);
}

/**
Expand Down Expand Up @@ -162,30 +202,64 @@ tensor3D SemiIsotropicBerendsenManostat::calculateMu() const
/**
* @brief calculate mu as scaling factor for Berendsen manostat (anisotropic)
*
* @details If a fixed axis is specified, that axis is not scaled (mu = 1.0)
* and the other axes are scaled independently
*
* @return tensor3D
*/
tensor3D AnisotropicBerendsenManostat::calculateMu() const
{
using enum FixedAxis;

const auto pxyz = diagonal(_pressureTensor);
const auto preFactor = _compressibility * _dt / _tau;

return diagonalMatrix(1.0 - preFactor * (_targetPressure - pxyz));
auto mu = 1.0 - preFactor * (_targetPressure - pxyz);

// 2D pressure coupling
if (_fixedAxis != NONE)
{
const auto fixedAxisIndex = static_cast<size_t>(_fixedAxis) - 1;
mu[fixedAxisIndex] = 1.0;
}

return diagonalMatrix(mu);
}

/**
* @brief calculate mu as scaling factor for Berendsen manostat (full
* anisotropic including angles)
*
* @details If a fixed axis is specified, the corresponding row and column
* are zeroed (no coupling with other axes) and the diagonal is set to 1.0
*
* @return tensor3D
*/
tensor3D FullAnisotropicBerendsenManostat::calculateMu() const
{
using enum FixedAxis;

const auto pTarget = diagonalMatrix(_targetPressure);
const auto preFactor = _compressibility * _dt / _tau;
const auto kronecker = kroneckerDeltaMatrix<double>();

auto mu = kronecker - preFactor * (pTarget - _pressureTensor);

// 2D full anisotropic: fix one axis and remove its coupling
if (_fixedAxis != NONE)
{
const auto fixedAxisIndex = static_cast<size_t>(_fixedAxis) - 1;

// Zero out the row and column of the fixed axis
for (size_t i = 0; i < 3; ++i)
{
mu[fixedAxisIndex][i] = 0.0;
mu[i][fixedAxisIndex] = 0.0;
}
// Set diagonal to 1.0 (no scaling)
mu[fixedAxisIndex][fixedAxisIndex] = 1.0;
}

rotateMu(mu);

return mu;
Expand Down
Loading
Loading