Skip to content

Commit

Permalink
Add ustruct volumetric penalty tests by testing rho and beta from g_v…
Browse files Browse the repository at this point in the history 
…ol_pen(). Also, add/update comments
  • Loading branch information
@aabrown100-git
aabrown100-git committed Jul 28, 2024
1 parent bfe86f9 commit 97131c0
Show file tree
Hide file tree
Showing 3 changed files with 238 additions and 18 deletions.
15 changes: 15 additions & 0 deletions Code/Source/svFSI/mat_models.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1434,6 +1434,21 @@ void get_tau(const ComMod& com_mod, const dmnType& lDmn, const double detF, cons
tauC = ctC * (he*c) * (rho0/detF);
}

/// @brief Compute rho, beta, drho/dp, dbeta/dp for volumetric penalty terms in
/// the ustruct formulation.
/// See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4
//
// ARGS:
// com_mod: ComMod object
// lDmn: dmnType object
// p: pressure
// ro: Solid density, rho
// bt: Isothermal compressibility coefficient, beta
// dro: Derivative of rho w.r.t. p
// dbt: Derivative of beta w.r.t. p
// Ja: Active strain Jacobian
//
// RETURNS: None, but updates ro, bt, dro, dbt
void g_vol_pen(const ComMod& com_mod, const dmnType& lDmn, const double p,
double& ro, double& bt, double& dro, double& dbt, const double Ja)
{
Expand Down
174 changes: 156 additions & 18 deletions tests/unitTests/test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -640,8 +640,20 @@ class QuadraticVolumetricPenaltyTest : public ::testing::Test {
}
};

// ------------------------------ STRUCT TESTS --------------------------------
// Test fixture class for STRUCT Quadratic penalty model
class STRUCT_QuadraticVolumetricPenaltyTest : public QuadraticVolumetricPenaltyTest {
protected:
void SetUp() override {
QuadraticVolumetricPenaltyTest::SetUp();

// Use struct
//TestQVP->ustruct = false;
}
};

// Test PK2 stress zero for F = I
TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressIdentityF) {
TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressIdentityF) {
//verbose = true; // Show values of S and S_ref

// Check identity F produces zero PK2 stress
Expand All @@ -653,7 +665,7 @@ TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressIdentityF) {
}

// Test PK2 stress zero for prescribed isochoric deformation
TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) {
TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) {
//verbose = true; // Show values of S and S_ref

// Check isochoric deformation produces zero PK2 stress
Expand All @@ -665,7 +677,7 @@ TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformati
}

// Test PK2 stress zero for random isochoric deformation
TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
//verbose = true; // Show values of S and S_ref

// Create random deformation gradient tensor
Expand All @@ -686,7 +698,7 @@ TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation)
}

// Test PK2 stress with finite difference for random F
TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressRandomF) {
TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressRandomF) {
//verbose = true; // Show values of S, dE, SdE and dPsi

// Compute reference PK2 stress with finite difference for random F
Expand All @@ -699,7 +711,7 @@ TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressRandomF) {
}

// Test PK2 stress consistent with strain energy for random F
TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
//verbose = true; // Show values of S, dE, SdE and dPsi

// Check random F produces consistent PK2 stress
Expand All @@ -708,14 +720,48 @@ TEST_F(QuadraticVolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
}

// Test material elasticity consistent with PK2 stress
TEST_F(QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) {
TEST_F(STRUCT_QuadraticVolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) {
//verbose = true; // Show values of CC, dE, CCdE and dS

// Check with random F
create_random_F(F);
TestQVP->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose);
}

// ------------------------------ USTRUCT TESTS --------------------------------
// Test fixture class for USTRUCT Quadratic penalty model
class USTRUCT_QuadraticVolumetricPenaltyTest : public QuadraticVolumetricPenaltyTest {
protected:
void SetUp() override {
QuadraticVolumetricPenaltyTest::SetUp();

// Use ustruct
//TestQVP->ustruct = true;
}
};

// Test rho, beta, drho/dp and dbeta/dp for random pressure
TEST_F(USTRUCT_QuadraticVolumetricPenaltyTest, TestRhoBeta) {
//verbose = true; // Show values of rho, beta, drho/dp and dbeta/dp and reference values

// Generate random pressure p between 100 and 1000
double p = getRandomDouble(100.0, 1000.0);

// Other parameters
double rho0 = 1000.0; // Reference density
double kappa = TestQVP->params.kappa; // Volumetric penalty parameter

// Compute reference values for rho, beta, drho/dp and dbeta/dp
// See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4
double rho_ref = rho0 / (1.0 - p/kappa);
double beta_ref = 1.0 / (kappa - p);
double drhodp_ref = -rho0 / pow(1.0 - p/kappa, 2) * (-1.0/kappa); // Derivative of rho with respect to p
double dbetadp_ref = -1.0 / pow(kappa - p, 2) * (-1.0); // Derivative of beta with respect to p

// Check rho, beta, drho/dp and dbeta/dp against reference values
TestQVP->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose);
}

// ----------------------------------------------------------------------------
// --------------------------- Simo-Taylor91 Volumetric Penalty Model ---------
// ----------------------------------------------------------------------------
Expand Down Expand Up @@ -753,8 +799,20 @@ class SimoTaylor91VolumetricPenaltyTest : public ::testing::Test {
}
};

// ------------------------------ STRUCT TESTS --------------------------------
// Test fixture class for STRUCT Simo-Taylor91 penalty model
class STRUCT_SimoTaylor91VolumetricPenaltyTest : public SimoTaylor91VolumetricPenaltyTest {
protected:
void SetUp() override {
SimoTaylor91VolumetricPenaltyTest::SetUp();

// Use struct
//TestST91->ustruct = false;
}
};

// Test PK2 stress zero for F = I
TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressIdentityF) {
TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressIdentityF) {
//verbose = true; // Show values of S and S_ref

// Check identity F produces zero PK2 stress
Expand All @@ -766,7 +824,7 @@ TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressIdentityF) {
}

// Test PK2 stress zero for prescribed isochoric deformation
TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) {
TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) {
//verbose = true; // Show values of S and S_ref

// Check isochoric deformation produces zero PK2 stress
Expand All @@ -778,7 +836,7 @@ TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeform
}

// Test PK2 stress zero for random isochoric deformation
TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
//verbose = true; // Show values of S and S_ref

// Create random deformation gradient tensor
Expand All @@ -799,7 +857,7 @@ TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformatio
}

// Test PK2 stress with finite difference for random F
TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomF) {
TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomF) {
//verbose = true; // Show values of S, dE, SdE and dPsi

// Compute reference PK2 stress with finite difference for random F
Expand All @@ -812,7 +870,7 @@ TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressRandomF) {
}

// Test PK2 stress consistent with strain energy for random F
TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
//verbose = true; // Show values of S, dE, SdE and dPsi

// Check random F produces consistent PK2 stress
Expand All @@ -821,14 +879,48 @@ TEST_F(SimoTaylor91VolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
}

// Test material elasticity consistent with PK2 stress
TEST_F(SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) {
TEST_F(STRUCT_SimoTaylor91VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) {
//verbose = true; // Show values of CC, dE, CCdE and dS

// Check with random F
create_random_F(F);
TestST91->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose);
}

// ------------------------------ USTRUCT TESTS --------------------------------
// Test fixture class for USTRUCT Simo-Taylor91 penalty model
class USTRUCT_SimoTaylor91VolumetricPenaltyTest : public SimoTaylor91VolumetricPenaltyTest {
protected:
void SetUp() override {
SimoTaylor91VolumetricPenaltyTest::SetUp();

// Use ustruct
//TestST91->ustruct = true;
}
};

// Test rho and beta values for random p
TEST_F(USTRUCT_SimoTaylor91VolumetricPenaltyTest, TestRhoBeta) {
//verbose = true; // Show values of rho, beta and rho_ref

// Generate random pressure p between 100 and 1000
double p = getRandomDouble(100.0, 1000.0);

// Other parameters
double rho0 = 1000.0; // Reference density
double kappa = TestST91->params.kappa; // Volumetric penalty parameter

// Compute reference values for rho, beta, drho/dp and dbeta/dp
// See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4
double rho_ref = rho0 / kappa * (sqrt(p*p + kappa*kappa) + p);
double beta_ref = 1.0 / sqrt(p*p + kappa*kappa);
double drhodp_ref = rho0 / kappa * (1.0/2.0 * (1.0/sqrt(p*p + kappa*kappa)) * 2.0*p + 1.0); // Derivative of rho with respect to p
double dbetadp_ref = -1.0/2.0 * 1.0/pow(p*p + kappa*kappa, 3.0/2.0) * 2.0*p; // Derivative of beta with respect to p

// Check rho, beta, drho/dp and dbeta/dp against reference values
TestST91->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose);
}

// ----------------------------------------------------------------------------
// --------------------------- Miehe94 Volumetric Penalty Model ---------------
// ----------------------------------------------------------------------------
Expand Down Expand Up @@ -865,8 +957,20 @@ class Miehe94VolumetricPenaltyTest : public ::testing::Test {
}
};

// ------------------------------ STRUCT TESTS --------------------------------
// Test fixture class for STRUCT Miehe94 penalty model
class STRUCT_Miehe94VolumetricPenaltyTest : public Miehe94VolumetricPenaltyTest {
protected:
void SetUp() override {
Miehe94VolumetricPenaltyTest::SetUp();

// Use struct
//TestM94->ustruct = false;
}
};

// Test PK2 stress zero for F = I
TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressIdentityF) {
TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressIdentityF) {
//verbose = true; // Show values of S and S_ref

// Check identity F produces zero PK2 stress
Expand All @@ -878,7 +982,7 @@ TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressIdentityF) {
}

// Test PK2 stress zero for prescribed isochoric deformation
TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) {
TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation) {
//verbose = true; // Show values of S and S_ref

// Check isochoric deformation produces zero PK2 stress
Expand All @@ -890,7 +994,7 @@ TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressPrescribedIsochoricDeformation
}

// Test PK2 stress zero for random isochoric deformation
TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
//verbose = true; // Show values of S and S_ref

// Create random deformation gradient tensor
Expand All @@ -911,7 +1015,7 @@ TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressRandomIsochoricDeformation) {
}

// Test PK2 stress with finite difference for random F
TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressRandomF) {
TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressRandomF) {
//verbose = true; // Show values of S, dE, SdE and dPsi

// Compute reference PK2 stress with finite difference for random F
Expand All @@ -924,7 +1028,7 @@ TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressRandomF) {
}

// Test PK2 stress consistent with strain energy for random F
TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
//verbose = true; // Show values of S, dE, SdE and dPsi

// Check random F produces consistent PK2 stress
Expand All @@ -933,14 +1037,48 @@ TEST_F(Miehe94VolumetricPenaltyTest, TestPK2StressConsistentRandomF) {
}

// Test material elasticity consistent with PK2 stress
TEST_F(Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) {
TEST_F(STRUCT_Miehe94VolumetricPenaltyTest, TestMaterialElasticityConsistentRandomF) {
//verbose = true; // Show values of CC, dE, CCdE and dS

// Check with random F
create_random_F(F);
TestM94->testMaterialElasticityConsistentWithPK2Stress(F, n_iter, rel_tol, abs_tol, delta, verbose);
}

// ------------------------------ USTRUCT TESTS --------------------------------
// Test fixture class for USTRUCT Miehe94 penalty model
class USTRUCT_Miehe94VolumetricPenaltyTest : public Miehe94VolumetricPenaltyTest {
protected:
void SetUp() override {
Miehe94VolumetricPenaltyTest::SetUp();

// Use ustruct
//TestM94->ustruct = true;
}
};

// Test rho and beta values for random p
TEST_F(USTRUCT_Miehe94VolumetricPenaltyTest, TestRhoBeta) {
//verbose = true; // Show values of rho, beta and rho_ref

// Generate random pressure p between 100 and 1000
double p = getRandomDouble(100.0, 1000.0);

// Other parameters
double rho0 = 1000.0; // Reference density
double kappa = TestM94->params.kappa; // Volumetric penalty parameter

// Compute reference values for rho, beta, drho/dp and dbeta/dp
// See ustruct paper (https://doi.org/10.1016/j.cma.2018.03.045) Section 2.4
double rho_ref = rho0 * (1.0 + p/kappa);
double beta_ref = 1.0 / (p + kappa);
double drhodp_ref = rho0 / kappa; // Derivative of rho with respect to p
double dbetadp_ref = -1.0 / pow(p + kappa, 2); // Derivative of beta with respect to p

// Check rho, beta, drho/dp and dbeta/dp against reference values
TestM94->testRhoBetaAgainstReference(p, rho0, rho_ref, beta_ref, drhodp_ref, dbetadp_ref, rel_tol, abs_tol, verbose);
}

// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------

Expand Down
Loading

0 comments on commit 97131c0

Please sign in to comment.