diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index aea6f047a..26b5218ec 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -47,6 +47,6 @@ jobs: export OMPI_ALLOW_RUN_AS_ROOT=1 export OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 - cd tests/automatic_tests + cd automatic_tests python3 run_automatic_tests.py diff --git a/.gitignore b/.gitignore index 04b65ca7d..8cd4badcb 100644 --- a/.gitignore +++ b/.gitignore @@ -63,6 +63,7 @@ main integratedFields.txt output.txt error.txt +CTestCostData.txt # Checkpoint/restart files restart.* diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/CMakeLists.txt b/automatic_tests/CHAC_anisotropyRegularized/CMakeLists.txt similarity index 94% rename from tests/automatic_tests/CHAC_anisotropyRegularized/CMakeLists.txt rename to automatic_tests/CHAC_anisotropyRegularized/CMakeLists.txt index 3b7e5c118..82b46ea0a 100644 --- a/tests/automatic_tests/CHAC_anisotropyRegularized/CMakeLists.txt +++ b/automatic_tests/CHAC_anisotropyRegularized/CMakeLists.txt @@ -136,8 +136,8 @@ if(EXISTS "nucleation.cc") endif() # Set location of files -include_directories(${CMAKE_SOURCE_DIR}/../../../include) -include_directories(${CMAKE_SOURCE_DIR}/../../../src) +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) include_directories(${CMAKE_SOURCE_DIR}) # Set the location of the main.cc file @@ -154,12 +154,12 @@ if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") add_executable(main_debug ${TARGET_SRC}) set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) deal_ii_setup_target(main_debug DEBUG) - target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-debug.a) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) endif() if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") add_executable(main_release ${TARGET_SRC}) set_property(TARGET main_release PROPERTY OUTPUT_NAME main) deal_ii_setup_target(main_release RELEASE) - target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-release.a) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) endif() \ No newline at end of file diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/ICs_and_BCs.cc b/automatic_tests/CHAC_anisotropyRegularized/ICs_and_BCs.cc similarity index 100% rename from tests/automatic_tests/CHAC_anisotropyRegularized/ICs_and_BCs.cc rename to automatic_tests/CHAC_anisotropyRegularized/ICs_and_BCs.cc diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/customPDE.h b/automatic_tests/CHAC_anisotropyRegularized/customPDE.h similarity index 100% rename from tests/automatic_tests/CHAC_anisotropyRegularized/customPDE.h rename to automatic_tests/CHAC_anisotropyRegularized/customPDE.h diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/equations.cc b/automatic_tests/CHAC_anisotropyRegularized/equations.cc similarity index 100% rename from tests/automatic_tests/CHAC_anisotropyRegularized/equations.cc rename to automatic_tests/CHAC_anisotropyRegularized/equations.cc diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/gold_integratedFields.txt b/automatic_tests/CHAC_anisotropyRegularized/gold_integratedFields.txt similarity index 100% rename from tests/automatic_tests/CHAC_anisotropyRegularized/gold_integratedFields.txt rename to automatic_tests/CHAC_anisotropyRegularized/gold_integratedFields.txt diff --git a/automatic_tests/CHAC_anisotropyRegularized/parameters.prm b/automatic_tests/CHAC_anisotropyRegularized/parameters.prm new file mode 100644 index 000000000..0e5d8a9e3 --- /dev/null +++ b/automatic_tests/CHAC_anisotropyRegularized/parameters.prm @@ -0,0 +1,77 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +set Domain size X = 100 +set Domain size Y = 100 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 1 +set Subdivisions Y = 1 +set Subdivisions Z = 1 + +set Refine factor = 7 + +set Element degree = 1 + +# ================================================================================= +# Set the adaptive mesh refinement parameters +# ================================================================================= +set Mesh adaptivity = true + +set Max refinement level = 7 +set Min refinement level = 4 + +set Steps between remeshing operations = 1000 + +subsection Refinement criterion: n + set Criterion type = VALUE + set Value lower bound = 0.001 + set Value upper bound = 0.999 +end + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 5.0e-2 + +set Number of time steps = 5000 + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable c = NATURAL +set Boundary condition for variable n = NATURAL +set Boundary condition for variable biharm = NATURAL + +# ================================================================================= +# Set the model constants +# ================================================================================= +# The CH mobility, McV in equations.h +set Model constant McV = 1.0, DOUBLE + +# The AC mobility, MnV in equations.h +set Model constant MnV = 0.1, DOUBLE + +# Anisotropy parameter +set Model constant epsilonM = 0.2, DOUBLE + +# Regularization parameter +set Model constant delta2 = 1.0, DOUBLE + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 10 + +set Skip print steps = 1000 diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/postprocess.cc b/automatic_tests/CHAC_anisotropyRegularized/postprocess.cc similarity index 100% rename from tests/automatic_tests/CHAC_anisotropyRegularized/postprocess.cc rename to automatic_tests/CHAC_anisotropyRegularized/postprocess.cc diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/CMakeLists.txt b/automatic_tests/CHiMaD_benchmark6a/CMakeLists.txt similarity index 94% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/CMakeLists.txt rename to automatic_tests/CHiMaD_benchmark6a/CMakeLists.txt index 3b7e5c118..82b46ea0a 100644 --- a/tests/automatic_tests/coupledCahnHilliardAllenCahn/CMakeLists.txt +++ b/automatic_tests/CHiMaD_benchmark6a/CMakeLists.txt @@ -136,8 +136,8 @@ if(EXISTS "nucleation.cc") endif() # Set location of files -include_directories(${CMAKE_SOURCE_DIR}/../../../include) -include_directories(${CMAKE_SOURCE_DIR}/../../../src) +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) include_directories(${CMAKE_SOURCE_DIR}) # Set the location of the main.cc file @@ -154,12 +154,12 @@ if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") add_executable(main_debug ${TARGET_SRC}) set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) deal_ii_setup_target(main_debug DEBUG) - target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-debug.a) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) endif() if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") add_executable(main_release ${TARGET_SRC}) set_property(TARGET main_release PROPERTY OUTPUT_NAME main) deal_ii_setup_target(main_release RELEASE) - target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-release.a) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) endif() \ No newline at end of file diff --git a/automatic_tests/CHiMaD_benchmark6a/ICs_and_BCs.cc b/automatic_tests/CHiMaD_benchmark6a/ICs_and_BCs.cc new file mode 100644 index 000000000..8617db09c --- /dev/null +++ b/automatic_tests/CHiMaD_benchmark6a/ICs_and_BCs.cc @@ -0,0 +1,81 @@ +// =========================================================================== +// FUNCTION FOR INITIAL CONDITIONS +// =========================================================================== + +template +void +customPDE::setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) +{ + // --------------------------------------------------------------------- + // ENTER THE INITIAL CONDITIONS HERE + // --------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the initial condition for each variable + // according to its variable index + + // The initial condition is a set of overlapping circles/spheres defined + // by a hyperbolic tangent function. The center of each circle/sphere is + // given by "center" and its radius is given by "radius". + + if (index == 0) + { + double x = p[0]; + double y = p[1]; + double c0 = 0.5; + double c1 = 0.04; + + double t1 = std::cos(0.2 * x) * std::cos(0.11 * y); + double t2 = std::cos(0.13 * x) * std::cos(0.087 * y) * std::cos(0.13 * x) * + std::cos(0.087 * y); + double t3 = std::cos(0.025 * x - 0.15 * y) * std::cos(0.07 * x - 0.02 * y); + + scalar_IC = c0 + c1 * (t1 + t2 + t3); + } + else + { + scalar_IC = 0.0; + } + + // ===================================================================== +} + +// =========================================================================== +// FUNCTION FOR NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS +// =========================================================================== + +template +void +customPDE::setNonUniformDirichletBCs( + [[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) +{ + // -------------------------------------------------------------------------- + // ENTER THE NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS HERE + // -------------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the boundary condition for each variable + // according to its variable index. This function can be left blank if there + // are no non-uniform Dirichlet boundary conditions. For BCs that change in + // time, you can access the current time through the variable "time". The + // boundary index can be accessed via the variable "direction", which starts + // at zero and uses the same order as the BC specification in parameters.in + // (i.e. left = 0, right = 1, bottom = 2, top = 3, front = 4, back = 5). + + if (index == 2) + { + if (direction == 1) + { + double y = p[1]; + scalar_BC = std::sin(y / 7.0); + } + } + + // ------------------------------------------------------------------------- +} diff --git a/automatic_tests/CHiMaD_benchmark6a/customPDE.h b/automatic_tests/CHiMaD_benchmark6a/customPDE.h new file mode 100644 index 000000000..5f4e54371 --- /dev/null +++ b/automatic_tests/CHiMaD_benchmark6a/customPDE.h @@ -0,0 +1,98 @@ +#include + +using namespace dealii; + +template +class customPDE : public MatrixFreePDE +{ +public: + // Constructor + customPDE(userInputParameters _userInputs) + : MatrixFreePDE(_userInputs) + , userInputs(_userInputs) {}; + + // Function to set the initial conditions (in ICs_and_BCs.h) + void + setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) override; + + // Function to set the non-uniform Dirichlet boundary conditions (in + // ICs_and_BCs.h) + void + setNonUniformDirichletBCs([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) override; + +private: +#include + + const userInputParameters userInputs; + + // Function to set the RHS of the governing equations for explicit time + // dependent equations (in equations.h) + void + explicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the RHS of the governing equations for all other equations + // (in equations.h) + void + nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the LHS of the governing equations (in equations.h) + void + equationLHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + +// Function to set postprocessing expressions (in postprocess.h) +#ifdef POSTPROCESS_FILE_EXISTS + void + postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; +#endif + +// Function to set the nucleation probability (in nucleation.h) +#ifdef NUCLEATION_FILE_EXISTS + double + getNucleationProbability([[maybe_unused]] variableValueContainer variable_value, + [[maybe_unused]] double dV) const override; +#endif + + // ================================================================ + // Methods specific to this subclass + // ================================================================ + + // ================================================================ + // Model constants specific to this subclass + // ================================================================ + + double McV = userInputs.get_model_constant_double("McV"); + double KcV = userInputs.get_model_constant_double("KcV"); + double rho = userInputs.get_model_constant_double("rho"); + double c_alpha = userInputs.get_model_constant_double("c_alpha"); + double c_beta = userInputs.get_model_constant_double("c_beta"); + double k = userInputs.get_model_constant_double("k"); + double epsilon = userInputs.get_model_constant_double("epsilon"); + + // ================================================================ +}; diff --git a/automatic_tests/CHiMaD_benchmark6a/equations.cc b/automatic_tests/CHiMaD_benchmark6a/equations.cc new file mode 100644 index 000000000..7e04322c0 --- /dev/null +++ b/automatic_tests/CHiMaD_benchmark6a/equations.cc @@ -0,0 +1,172 @@ +// ================================================================================= +// Set the attributes of the primary field variables +// ================================================================================= +// This function sets attributes for each variable/equation in the app. The +// attributes are set via standardized function calls. The first parameter for +// each function call is the variable index (starting at zero). The first set of +// variable/equation attributes are the variable name (any string), the variable +// type (SCALAR/VECTOR), and the equation type (EXPLICIT_TIME_DEPENDENT/ +// TIME_INDEPENDENT/AUXILIARY). The next set of attributes describe the +// dependencies for the governing equation on the values and derivatives of the +// other variables for the value term and gradient term of the RHS and the LHS. +// The final pair of attributes determine whether a variable represents a field +// that can nucleate and whether the value of the field is needed for nucleation +// rate calculations. + +void +variableAttributeLoader::loadVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "c"); + set_variable_type(0, SCALAR); + set_variable_equation_type(0, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(0, "c"); + set_dependencies_gradient_term_RHS(0, "grad(mu)"); + + // Variable 1 + set_variable_name(1, "mu"); + set_variable_type(1, SCALAR); + set_variable_equation_type(1, AUXILIARY); + + set_dependencies_value_term_RHS(1, "c, phi"); + set_dependencies_gradient_term_RHS(1, "grad(c)"); + + // Variable 2 + set_variable_name(2, "phi"); + set_variable_type(2, SCALAR); + set_variable_equation_type(2, TIME_INDEPENDENT); + + set_dependencies_value_term_RHS(2, "c"); + set_dependencies_gradient_term_RHS(2, "grad(phi)"); + set_dependencies_value_term_LHS(2, ""); + set_dependencies_gradient_term_LHS(2, "grad(change(phi))"); +} + +// ============================================================================================= +// explicitEquationRHS (needed only if one or more equation is explict time +// dependent) +// ============================================================================================= +// This function calculates the right-hand-side of the explicit time-dependent +// equations for each variable. It takes "variable_list" as an input, which is a +// list of the value and derivatives of each of the variables at a specific +// quadrature point. The (x,y,z) location of that quadrature point is given by +// "q_point_loc". The function outputs two terms to variable_list -- one +// proportional to the test function and one proportional to the gradient of the +// test function. The index for each variable in this list corresponds to the +// index given at the top of this file. + +template +void +customPDE::explicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The concentration and its derivatives + scalarvalueType c = variable_list.get_scalar_value(0); + + // The chemical potential and its derivatives + scalargradType mux = variable_list.get_scalar_gradient(1); + + // --- Setting the expressions for the terms in the governing equations --- + + // The terms in the equations + scalarvalueType eq_c = c; + scalargradType eqx_c = constV(-McV * userInputs.dtValue) * mux; + + // --- Submitting the terms for the governing equations --- + + variable_list.set_scalar_value_term_RHS(0, eq_c); + variable_list.set_scalar_gradient_term_RHS(0, eqx_c); +} + +// ============================================================================================= +// nonExplicitEquationRHS (needed only if one or more equation is time +// independent or auxiliary) +// ============================================================================================= +// This function calculates the right-hand-side of all of the equations that are +// not explicit time-dependent equations. It takes "variable_list" as an input, +// which is a list of the value and derivatives of each of the variables at a +// specific quadrature point. The (x,y,z) location of that quadrature point is +// given by "q_point_loc". The function outputs two terms to variable_list -- +// one proportional to the test function and one proportional to the gradient of +// the test function. The index for each variable in this list corresponds to +// the index given at the top of this file. + +template +void +customPDE::nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + scalarvalueType c = variable_list.get_scalar_value(0); + scalargradType cx = variable_list.get_scalar_gradient(0); + + // The electric potential and its derivatives + scalarvalueType phi = variable_list.get_scalar_value(2); + scalargradType phix = variable_list.get_scalar_gradient(2); + + // --- Setting the expressions for the terms in the governing equations --- + + // The derivative of the local chemical free energy + scalarvalueType fcV = + 2.0 * rho * (c - c_alpha) * (c_beta - c) * (c_alpha + c_beta - 2.0 * c); + + // The derivative of the local electrostatic free energy + scalarvalueType fphiV = k * phi; + + scalarvalueType eq_mu = fcV + fphiV; + scalargradType eqx_mu = constV(KcV) * cx; + scalarvalueType eq_phi = -constV(-k / epsilon) * c; + scalargradType eqx_phi = -phix; + + // --- Submitting the terms for the governing equations --- + + variable_list.set_scalar_value_term_RHS(1, eq_mu); + variable_list.set_scalar_gradient_term_RHS(1, eqx_mu); + + variable_list.set_scalar_value_term_RHS(2, eq_phi); + variable_list.set_scalar_gradient_term_RHS(2, eqx_phi); +} + +// ============================================================================================= +// equationLHS (needed only if at least one equation is time independent) +// ============================================================================================= +// This function calculates the left-hand-side of time-independent equations. It +// takes "variable_list" as an input, which is a list of the value and +// derivatives of each of the variables at a specific quadrature point. The +// (x,y,z) location of that quadrature point is given by "q_point_loc". The +// function outputs two terms to variable_list -- one proportional to the test +// function and one proportional to the gradient of the test function -- for the +// left-hand-side of the equation. The index for each variable in this list +// corresponds to the index given at the top of this file. If there are multiple +// elliptic equations, conditional statements should be sed to ensure that the +// correct residual is being submitted. The index of the field being solved can +// be accessed by "this->currentFieldIndex". + +template +void +customPDE::equationLHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // grad(delta phi) + scalargradType Dphix = variable_list.get_change_in_scalar_gradient(2); + + // --- Setting the expressions for the terms in the governing equations --- + + scalargradType eqx_Dphi = Dphix; + + // --- Submitting the terms for the governing equations --- + + variable_list.set_scalar_gradient_term_LHS(2, eqx_Dphi); +} diff --git a/automatic_tests/CHiMaD_benchmark6a/gold_integratedFields.txt b/automatic_tests/CHiMaD_benchmark6a/gold_integratedFields.txt new file mode 100644 index 000000000..aa1851e44 --- /dev/null +++ b/automatic_tests/CHiMaD_benchmark6a/gold_integratedFields.txt @@ -0,0 +1,11 @@ +0 f_tot 183.0435848 +0.4 f_tot 182.6412281 +0.8 f_tot 182.1957345 +1.2 f_tot 181.6864768 +1.6 f_tot 181.1100944 +2 f_tot 180.455664 +2.4 f_tot 179.7174222 +2.8 f_tot 178.8862819 +3.2 f_tot 177.958548 +3.6 f_tot 176.9306608 +4 f_tot 175.8053599 diff --git a/automatic_tests/CHiMaD_benchmark6a/parameters.prm b/automatic_tests/CHiMaD_benchmark6a/parameters.prm new file mode 100644 index 000000000..1d61e18cb --- /dev/null +++ b/automatic_tests/CHiMaD_benchmark6a/parameters.prm @@ -0,0 +1,92 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +set Domain size X = 100 +set Domain size Y = 100 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 1 +set Subdivisions Y = 1 +set Subdivisions Z = 1 + +set Refine factor = 6 + +set Element degree = 2 + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 2.0e-4 + +set Number of time steps = 20000 + +# ================================================================================= +# Set the linear solver parameters +# ================================================================================= + +subsection Linear solver parameters: phi + set Tolerance type = ABSOLUTE_RESIDUAL + + set Tolerance value = 1e-5 + + set Maximum linear solver iterations = 10000 +end + +# ================================================================================= +# Set the nonlinear solver parameters +# ================================================================================= + +set Maximum nonlinear solver iterations = 100 + +subsection Nonlinear solver parameters: mu + set Tolerance type = ABSOLUTE_SOLUTION_CHANGE + set Tolerance value = 1e-5 + set Use backtracking line search damping = false + set Constant damping value = 1.0 + set Use Laplace's equation to determine the initial guess = false +end + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable c = NATURAL +set Boundary condition for variable mu = NATURAL +set Boundary condition for variable phi = DIRICHLET: 0.0, NON_UNIFORM_DIRICHLET, NATURAL, NATURAL, NATURAL, NATURAL + +# ================================================================================= +# Set the model constants +# ================================================================================= +# The mobility, McV in equations.h +set Model constant McV = 5.0, DOUBLE + +# The gradient energy coefficient, KcV in equations.h +set Model constant KcV = 2.0, DOUBLE + +# Parameters for the chemical and electrostatic components of the free energy + +#chemical +set Model constant rho = 5.0, DOUBLE +set Model constant c_alpha = 0.3, DOUBLE +set Model constant c_beta = 0.7, DOUBLE + +# electrostatic +set Model constant k = 0.09, DOUBLE +set Model constant epsilon = 90.0, DOUBLE + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 10 + +set Skip print steps = 1000 diff --git a/automatic_tests/CHiMaD_benchmark6a/postprocess.cc b/automatic_tests/CHiMaD_benchmark6a/postprocess.cc new file mode 100644 index 000000000..7292e9eb6 --- /dev/null +++ b/automatic_tests/CHiMaD_benchmark6a/postprocess.cc @@ -0,0 +1,83 @@ +// ============================================================================================= +// loadPostProcessorVariableAttributes: Set the attributes of the postprocessing +// variables +// ============================================================================================= +// This function is analogous to 'loadVariableAttributes' in 'equations.h', but +// for the postprocessing expressions. It sets the attributes for each +// postprocessing expression, including its name, whether it is a vector or +// scalar (only scalars are supported at present), its dependencies on other +// variables and their derivatives, and whether to calculate an integral of the +// postprocessed quantity over the entire domain. Note: this function is not a +// member of customPDE. + +void +variableAttributeLoader::loadPostProcessorVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "f_tot"); + set_variable_type(0, SCALAR); + + set_dependencies_value_term_RHS(0, "c, grad(c), phi"); + set_dependencies_gradient_term_RHS(0, ""); + + set_output_integral(0, true); +} + +// ============================================================================================= +// postProcessedFields: Set the postprocessing expressions +// ============================================================================================= +// This function is analogous to 'explicitEquationRHS' and +// 'nonExplicitEquationRHS' in equations.h. It takes in "variable_list" and +// "q_point_loc" as inputs and outputs two terms in the expression for the +// postprocessing variable -- one proportional to the test function and one +// proportional to the gradient of the test function. The index for each +// variable in this list corresponds to the index given at the top of this file +// (for submitting the terms) and the index in 'equations.h' for assigning the +// values/derivatives of the primary variables. + +template +void +customPDE::postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The concentration and its derivatives + scalarvalueType c = variable_list.get_scalar_value(0); + scalargradType cx = variable_list.get_scalar_gradient(0); + scalarvalueType phi = variable_list.get_scalar_value(2); + + // --- Setting the expressions for the terms in the postprocessing expressions + // --- + + scalarvalueType f_tot = constV(0.0); + + // The homogenous chemical free energy + scalarvalueType f_chem = + rho * (c - c_alpha) * (c - c_alpha) * (c_beta - c) * (c_beta - c); + + // The homogenous electrostaric free energy + scalarvalueType f_elec = 0.5 * k * c * phi; + + // The gradient free energy + scalarvalueType f_grad = constV(0.0); + + for (int i = 0; i < dim; i++) + { + for (int j = 0; j < dim; j++) + { + f_grad += constV(0.5 * KcV) * cx[i] * cx[j]; + } + } + + // The total free energy + f_tot = f_chem + f_elec + f_grad; + + // --- Submitting the terms for the postprocessing expressions --- + pp_variable_list.set_scalar_value_term_RHS(0, f_tot); +} diff --git a/tests/automatic_tests/allenCahn/CMakeLists.txt b/automatic_tests/allenCahn/CMakeLists.txt similarity index 94% rename from tests/automatic_tests/allenCahn/CMakeLists.txt rename to automatic_tests/allenCahn/CMakeLists.txt index 3b7e5c118..82b46ea0a 100644 --- a/tests/automatic_tests/allenCahn/CMakeLists.txt +++ b/automatic_tests/allenCahn/CMakeLists.txt @@ -136,8 +136,8 @@ if(EXISTS "nucleation.cc") endif() # Set location of files -include_directories(${CMAKE_SOURCE_DIR}/../../../include) -include_directories(${CMAKE_SOURCE_DIR}/../../../src) +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) include_directories(${CMAKE_SOURCE_DIR}) # Set the location of the main.cc file @@ -154,12 +154,12 @@ if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") add_executable(main_debug ${TARGET_SRC}) set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) deal_ii_setup_target(main_debug DEBUG) - target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-debug.a) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) endif() if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") add_executable(main_release ${TARGET_SRC}) set_property(TARGET main_release PROPERTY OUTPUT_NAME main) deal_ii_setup_target(main_release RELEASE) - target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-release.a) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) endif() \ No newline at end of file diff --git a/tests/automatic_tests/allenCahn/ICs_and_BCs.cc b/automatic_tests/allenCahn/ICs_and_BCs.cc similarity index 100% rename from tests/automatic_tests/allenCahn/ICs_and_BCs.cc rename to automatic_tests/allenCahn/ICs_and_BCs.cc diff --git a/tests/automatic_tests/allenCahn/customPDE.h b/automatic_tests/allenCahn/customPDE.h similarity index 100% rename from tests/automatic_tests/allenCahn/customPDE.h rename to automatic_tests/allenCahn/customPDE.h diff --git a/tests/automatic_tests/allenCahn/equations.cc b/automatic_tests/allenCahn/equations.cc similarity index 100% rename from tests/automatic_tests/allenCahn/equations.cc rename to automatic_tests/allenCahn/equations.cc diff --git a/tests/automatic_tests/allenCahn/gold_integratedFields.txt b/automatic_tests/allenCahn/gold_integratedFields.txt similarity index 100% rename from tests/automatic_tests/allenCahn/gold_integratedFields.txt rename to automatic_tests/allenCahn/gold_integratedFields.txt diff --git a/automatic_tests/allenCahn/parameters.prm b/automatic_tests/allenCahn/parameters.prm new file mode 100644 index 000000000..3e68a14c8 --- /dev/null +++ b/automatic_tests/allenCahn/parameters.prm @@ -0,0 +1,55 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +set Domain size X = 100 +set Domain size Y = 100 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 1 +set Subdivisions Y = 1 +set Subdivisions Z = 1 + +set Refine factor = 8 + +set Element degree = 1 + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 1.0e-2 + +set Number of time steps = 5000 + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 5 + +set Print timing information with output = true + +set Skip print steps = 1000 + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable n = NATURAL + +# ================================================================================= +# Set the model constants +# ================================================================================= +# The mobility, MnV in equations.cc +set Model constant MnV = 1.0, DOUBLE + +# The gradient energy coefficient, KnV in equations.cc +set Model constant KnV = 2.0, DOUBLE diff --git a/tests/automatic_tests/allenCahn/postprocess.cc b/automatic_tests/allenCahn/postprocess.cc similarity index 100% rename from tests/automatic_tests/allenCahn/postprocess.cc rename to automatic_tests/allenCahn/postprocess.cc diff --git a/tests/automatic_tests/application_debug_test.py b/automatic_tests/application_debug_test.py similarity index 96% rename from tests/automatic_tests/application_debug_test.py rename to automatic_tests/application_debug_test.py index 8f8b71063..586f0c0df 100644 --- a/tests/automatic_tests/application_debug_test.py +++ b/automatic_tests/application_debug_test.py @@ -36,12 +36,10 @@ def get_application_path(app_name): pwd = os.getcwd() # Check that we're in the automatic test directory - assert ( - "tests/automatic_tests" in pwd - ), "Current directory is not within 'tests/automatic_tests'" + assert "automatic_tests" in pwd, "Current directory is not within 'automatic_tests'" # Application path assuming file structure matches GitHub repo - app_path = pwd.replace("tests/automatic_tests", f"applications/{app_name}") + app_path = pwd.replace("automatic_tests", f"applications/{app_name}") return app_path @@ -101,7 +99,7 @@ def compile_and_run(app_name, new_parameter_file, test_dir): is the application name and the second is success or failure """ try: - # Navigate to tests/automatic_tests + # Navigate to automatic_tests os.chdir(test_dir) # Grab application path diff --git a/tests/automatic_tests/cahnHilliard/CMakeLists.txt b/automatic_tests/cahnHilliard/CMakeLists.txt similarity index 94% rename from tests/automatic_tests/cahnHilliard/CMakeLists.txt rename to automatic_tests/cahnHilliard/CMakeLists.txt index 3b7e5c118..82b46ea0a 100644 --- a/tests/automatic_tests/cahnHilliard/CMakeLists.txt +++ b/automatic_tests/cahnHilliard/CMakeLists.txt @@ -136,8 +136,8 @@ if(EXISTS "nucleation.cc") endif() # Set location of files -include_directories(${CMAKE_SOURCE_DIR}/../../../include) -include_directories(${CMAKE_SOURCE_DIR}/../../../src) +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) include_directories(${CMAKE_SOURCE_DIR}) # Set the location of the main.cc file @@ -154,12 +154,12 @@ if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") add_executable(main_debug ${TARGET_SRC}) set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) deal_ii_setup_target(main_debug DEBUG) - target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-debug.a) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) endif() if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") add_executable(main_release ${TARGET_SRC}) set_property(TARGET main_release PROPERTY OUTPUT_NAME main) deal_ii_setup_target(main_release RELEASE) - target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-release.a) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) endif() \ No newline at end of file diff --git a/tests/automatic_tests/cahnHilliard/ICs_and_BCs.cc b/automatic_tests/cahnHilliard/ICs_and_BCs.cc similarity index 100% rename from tests/automatic_tests/cahnHilliard/ICs_and_BCs.cc rename to automatic_tests/cahnHilliard/ICs_and_BCs.cc diff --git a/tests/automatic_tests/cahnHilliard/customPDE.h b/automatic_tests/cahnHilliard/customPDE.h similarity index 100% rename from tests/automatic_tests/cahnHilliard/customPDE.h rename to automatic_tests/cahnHilliard/customPDE.h diff --git a/tests/automatic_tests/cahnHilliard/equations.cc b/automatic_tests/cahnHilliard/equations.cc similarity index 100% rename from tests/automatic_tests/cahnHilliard/equations.cc rename to automatic_tests/cahnHilliard/equations.cc diff --git a/tests/automatic_tests/cahnHilliard/gold_integratedFields.txt b/automatic_tests/cahnHilliard/gold_integratedFields.txt similarity index 100% rename from tests/automatic_tests/cahnHilliard/gold_integratedFields.txt rename to automatic_tests/cahnHilliard/gold_integratedFields.txt diff --git a/automatic_tests/cahnHilliard/parameters.prm b/automatic_tests/cahnHilliard/parameters.prm new file mode 100644 index 000000000..69541cd3d --- /dev/null +++ b/automatic_tests/cahnHilliard/parameters.prm @@ -0,0 +1,53 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +set Domain size X = 100 +set Domain size Y = 100 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 1 +set Subdivisions Y = 1 +set Subdivisions Z = 1 + +set Refine factor = 6 + +set Element degree = 2 + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 1.0e-3 + +set Number of time steps = 5000 + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable c = NATURAL +set Boundary condition for variable mu = NATURAL + +# ================================================================================= +# Set the model constants +# The mobility, McV in equations.h +set Model constant McV = 1.0, DOUBLE + +# The gradient energy coefficient, KcV in equations.h +set Model constant KcV = 1.5, DOUBLE + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 10 + +set Skip print steps = 1000 diff --git a/tests/automatic_tests/cahnHilliard/postprocess.cc b/automatic_tests/cahnHilliard/postprocess.cc similarity index 100% rename from tests/automatic_tests/cahnHilliard/postprocess.cc rename to automatic_tests/cahnHilliard/postprocess.cc diff --git a/automatic_tests/corrosion_microgalvanic/CMakeLists.txt b/automatic_tests/corrosion_microgalvanic/CMakeLists.txt new file mode 100644 index 000000000..82b46ea0a --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/CMakeLists.txt @@ -0,0 +1,165 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/automatic_tests/corrosion_microgalvanic/ICs_and_BCs.cc b/automatic_tests/corrosion_microgalvanic/ICs_and_BCs.cc new file mode 100644 index 000000000..63501011d --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/ICs_and_BCs.cc @@ -0,0 +1,124 @@ +// =========================================================================== +// FUNCTION FOR INITIAL CONDITIONS +// =========================================================================== +#include +#include +using namespace std; + +template +void +customPDE::setInitialCondition( + [[maybe_unused]] const dealii::Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] dealii::Vector &vector_IC) +{ + // --------------------------------------------------------------------- + // ENTER THE INITIAL CONDITIONS HERE FOR SCALAR FIELDS + // --------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the initial condition for each variable + // according to its variable index + + // The initial condition is a set of overlapping circles/spheres defined + // by a hyperbolic tangent function. The center of each circle/sphere is + // given by "center" and its radius is given by "radius". + + /* + List of scalar variables + Index: 0; var. nAnodic + Index: 1; var. muAnodic + Index: 2; var. nCathodic + Index: 3; var. muCathodic + Index: 4; var. psi + Index: 5; var. mupsi + Index: 6; var. Phi + Index: 7; var. irxn + */ + + // xCoordinate of the point considered + double posx = p[0]; + + // yCoordinate of the point considered + double posy = p[1]; + + // zCoordinate of the point considered + // double posz = p[2]; + + const double SYSTEM_WIDTH = userInputs.domain_size[0]; + const double SYSTEM_HEIGHT = userInputs.domain_size[1]; + + double hyperTan = tanh((abs(posx - 1 * SYSTEM_WIDTH) - 1 * cathodeThickness) / deltaV); + double height = 0.5 * SYSTEM_HEIGHT; + + double orderParaCathodic = 0.5 * (1 - hyperTan); + double orderParaAnodic = 0.5 * (1 + hyperTan); + double domainParaPsi = 0.5 * (1 + tanh((posy - height) / deltaV)); + + orderParaCathodic = orderParaCathodic * (1 - domainParaPsi); + orderParaAnodic = orderParaAnodic * (1 - domainParaPsi); + + if (index == 0) + { // order parameter for the anodic phase, nAnodic + + scalar_IC = orderParaAnodic; + } + else if (index == 2) + { // order parameter for the cathodic phase, nCathodic + + scalar_IC = orderParaCathodic; + } + else if (index == 4) + { // electrolyte domain parameter, psi + + scalar_IC = domainParaPsi; + } + else if (index == 6) + { // electrolyte potential Phi + + scalar_IC = initialGuessForPhi; + } + else + { // all other variable fields + + scalar_IC = 0.0; + } + // --------------------------------------------------------------------- +} + +// =========================================================================== +// FUNCTION FOR NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS +// =========================================================================== + +template +void +customPDE::setNonUniformDirichletBCs( + [[maybe_unused]] const dealii::Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] dealii::Vector &vector_BC) +{ + // -------------------------------------------------------------------------- + // ENTER THE NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS HERE + // -------------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the boundary condition for each variable + // according to its variable index. This function can be left blank if there + // are no non-uniform Dirichlet boundary conditions. For BCs that change in + // time, you can access the current time through the variable "time". The + // boundary index can be accessed via the variable "direction", which starts + // at zero and uses the same order as the BC specification in parameters.in + // (i.e. left = 0, right = 1, bottom = 2, top = 3, front = 4, back = 5). + /* +if (index == 2){ + if (direction == 1){ + double x=p[0]; + double y=p[1]; + scalar_BC=sin(y/7.0); + } +} + */ + // ------------------------------------------------------------------------- +} diff --git a/automatic_tests/corrosion_microgalvanic/customPDE.h b/automatic_tests/corrosion_microgalvanic/customPDE.h new file mode 100644 index 000000000..a4edcbf5d --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/customPDE.h @@ -0,0 +1,112 @@ +#include +#include +#include +using namespace std; + +template +class customPDE : public MatrixFreePDE +{ +public: + // Constructor + customPDE(userInputParameters _userInputs) + : MatrixFreePDE(_userInputs) + , userInputs(_userInputs) {}; + + // Function to set the initial conditions (in ICs_and_BCs.h) + void + setInitialCondition([[maybe_unused]] const dealii::Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] dealii::Vector &vector_IC) override; + + // Function to set the non-uniform Dirichlet + // boundary conditions (in ICs_and_BCs.h) + void + setNonUniformDirichletBCs([[maybe_unused]] const dealii::Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] dealii::Vector &vector_BC) override; + +private: +#include + + const userInputParameters userInputs; + + // Function to set the RHS of the governing equations + // for explicit time dependent equations (in equations.h) + + void + explicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the RHS of the governing equations + // for all other equations (in equations.h) + void + nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the LHS of the governing equations (in equations.h) + void + equationLHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set postprocessing expressions (in postprocess.h) + +#ifdef POSTPROCESS_FILE_EXISTS + void + postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; +#endif + +// Function to set the nucleation probability (in nucleation.h) +// Not useful for this application +#ifdef NUCLEATION_FILE_EXISTS + double + getNucleationProbability([[maybe_unused]] variableValueContainer variable_value, + [[maybe_unused]] double dV) const override; +#endif + + // ================================================================ + // Model constants specific to this subclass + // ================================================================ + + double VMV = userInputs.get_model_constant_double("VMV"); + double zMV = userInputs.get_model_constant_double("zMV"); + double epssqV = userInputs.get_model_constant_double("epssqV"); + double EcorrAnodic = userInputs.get_model_constant_double("EcorrAnodic"); + double EcorrCathodic = userInputs.get_model_constant_double("EcorrCathodic"); + double AAnodic = userInputs.get_model_constant_double("AAnodic"); + double ACathodic = userInputs.get_model_constant_double("ACathodic"); + double VsV = userInputs.get_model_constant_double("VsV"); + double lthresh = userInputs.get_model_constant_double("lthresh"); + double gamma = userInputs.get_model_constant_double("gamma"); + double kappa = userInputs.get_model_constant_double("kappa"); + double i0Anodic = userInputs.get_model_constant_double("i0Anodic"); + double i0Cathodic = userInputs.get_model_constant_double("i0Cathodic"); + double iMax = userInputs.get_model_constant_double("iMax"); + double cathodeThickness = userInputs.get_model_constant_double("cathodeThickness"); + double tStepStartForV = userInputs.get_model_constant_double("tStepStartForV"); + double initialGuessForPhi = userInputs.get_model_constant_double("guessValPhi"); + double FarC = 96485.33289; // Faraday constant [C/mol] + double RV = 8.314; // Gas constant [J/(mol K)] + double deltaV = sqrt(2.0 * epssqV); // Equilibrium interface half-width + double MconstV = 2.0 * (VMV / (zMV * FarC)) * deltaV; // Constant prefactor for mobility + + // ================================================================ +}; diff --git a/automatic_tests/corrosion_microgalvanic/equations.cc b/automatic_tests/corrosion_microgalvanic/equations.cc new file mode 100644 index 000000000..b0b827fd3 --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/equations.cc @@ -0,0 +1,470 @@ +// ================================================================================= +// Set the attributes of the primary field variables +// ================================================================================= +// This function sets attributes for each variable/equation in the app. The +// attributes are set via standardized function calls. The first parameter for +// each function call is the variable index (starting at zero). The first set of +// variable/equation attributes are the variable name (any string), the variable +// type (SCALAR/VECTOR), and the equation type (EXPLICIT_TIME_DEPENDENT/ +// TIME_INDEPENDENT/AUXILIARY). The next set of attributes describe the +// dependencies for the governing equation on the values and derivatives of the +// other variables for the value term and gradient term of the RHS and the LHS. +// The final pair of attributes determine whether a variable represents a field +// that can nucleate and whether the value of the field is needed for nucleation +// rate calculations. +#include +#include +using namespace std; + +void +variableAttributeLoader::loadVariableAttributes() +{ + // Declaring all the model variables as follows: + // domain parameter for the electrolyte - psi + // order parameters for the anodic phase, nAnodic, and the cathodic phase, + // nCathodic + // chemical potentials for the anodic phase, muAnodic, + // and the cathodic phase, + // muCathodic + // electrolyte potential - Phi + // reaction current density - irxn + // interpolation factor for the anodic phase, xiAnodic + + // Variable 0 + set_variable_name(0, "nAnodic"); + set_variable_type(0, SCALAR); + set_variable_equation_type(0, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(0, "nAnodic, irxn, grad(psi), Phi, xiAnodic"); + set_dependencies_gradient_term_RHS(0, "nAnodic, grad(muAnodic), irxn, Phi, xiAnodic"); + + // Variable 1 + set_variable_name(1, "muAnodic"); + set_variable_type(1, SCALAR); + set_variable_equation_type(1, AUXILIARY); + + set_dependencies_value_term_RHS(1, "nAnodic, nCathodic, psi"); + set_dependencies_gradient_term_RHS(1, "grad(nAnodic)"); + + // Variable 2 + set_variable_name(2, "nCathodic"); + set_variable_type(2, SCALAR); + set_variable_equation_type(2, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(2, "nCathodic"); + set_dependencies_gradient_term_RHS(2, "nCathodic, grad(muCathodic), irxn, Phi"); + + // Variable 3 + set_variable_name(3, "muCathodic"); + set_variable_type(3, SCALAR); + set_variable_equation_type(3, AUXILIARY); + + set_dependencies_value_term_RHS(3, "nCathodic, nAnodic, psi"); + set_dependencies_gradient_term_RHS(3, "grad(nCathodic)"); + + // Variable 4 + set_variable_name(4, "psi"); + set_variable_type(4, SCALAR); + set_variable_equation_type(4, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(4, "psi, irxn, grad(psi), Phi, xiAnodic"); + set_dependencies_gradient_term_RHS(4, "psi, grad(mupsi), irxn, Phi, xiAnodic"); + + // Variable 5 + set_variable_name(5, "mupsi"); + set_variable_type(5, SCALAR); + set_variable_equation_type(5, AUXILIARY); + + set_dependencies_value_term_RHS(5, "nAnodic, nCathodic, psi"); + set_dependencies_gradient_term_RHS(5, "grad(psi)"); + + // Variable 6 + set_variable_name(6, "Phi"); + set_variable_type(6, SCALAR); + set_variable_equation_type(6, TIME_INDEPENDENT); + + set_dependencies_value_term_LHS( + 6, + "nAnodic, nCathodic, grad(psi), change(Phi), Phi, xiAnodic"); + set_dependencies_gradient_term_LHS(6, "psi, grad(change(Phi))"); + set_dependencies_value_term_RHS(6, "grad(psi), irxn, xiAnodic"); + set_dependencies_gradient_term_RHS(6, "psi, grad(Phi)"); + + // Variable 7 + set_variable_name(7, "irxn"); + set_variable_type(7, SCALAR); + set_variable_equation_type(7, AUXILIARY); + + set_dependencies_value_term_RHS(7, "nCathodic, nAnodic, Phi, xiAnodic"); + set_dependencies_gradient_term_RHS(7, ""); + + // Variable 8 + set_variable_name(8, "xiAnodic"); + set_variable_type(8, SCALAR); + set_variable_equation_type(8, AUXILIARY); + + set_dependencies_value_term_RHS(8, "nAnodic, nCathodic"); +} + +// ============================================================================================= +// explicitEquationRHS (needed only if one or more equation is explict time +// dependent) +// ============================================================================================= +// This function calculates the right-hand-side of the explicit time-dependent +// equations for each variable. It takes "variable_list" as an input, which is a +// list of the value and derivatives of each of the variables at a specific +// quadrature point. The (x,y,z) location of that quadrature point is given by +// "q_point_loc". The function outputs two terms to variable_list -- one +// proportional to the test function and one proportional to the gradient of the +// test function. The index for each variable in this list corresponds to the +// index given at the top of this file. + +template +void +customPDE::explicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Parameters in the explicit equations can be set here --- + + // Timestep + scalarvalueType delt = constV(userInputs.dtValue); + + // The order parameter of the anodic phase + scalarvalueType nAnodic = variable_list.get_scalar_value(0); + + // The gradient of the chemical potential of the anodic phase + scalargradType muAnodicx = variable_list.get_scalar_gradient(1); + + // The order parameter of the cathodic phase + scalarvalueType nCathodic = variable_list.get_scalar_value(2); + + // The gradient of the chemical potential of the cathodic phase + scalargradType muCathodicx = variable_list.get_scalar_gradient(3); + + // The domain parameter of the electrolyte phase and its derivatives + scalarvalueType psi = variable_list.get_scalar_value(4); + scalargradType psix = variable_list.get_scalar_gradient(4); + + // The derivative of the chemical potential of the electrolyte phase + scalargradType mupsix = variable_list.get_scalar_gradient(5); + + // The electrolyte potential + scalarvalueType Phi = variable_list.get_scalar_value(6); + + // The interpolation factor for the anodic phase + scalarvalueType xiAnodic = variable_list.get_scalar_value(8); + + // Ensuring that the domain and order parameters remain within [0 1] + + // Electrolyte + psi = min(psi, constV(1.0)); + psi = max(psi, constV(0.0)); + + // Anodic phase + nAnodic = min(nAnodic, constV(1.0)); + nAnodic = max(nAnodic, constV(0.0)); + + // Cathodic phase + nCathodic = min(nCathodic, constV(1.0)); + nCathodic = max(nCathodic, constV(0.0)); + + // --- Calculation of terms needed in multiple expressions --- + + // Magnifude of the derivative of psi + scalarvalueType magpsix = constV(0.0); + for (int i = 0; i < dim; i++) + { + magpsix += psix[i] * psix[i]; + } + magpsix = sqrt(magpsix); + + // --- Calculation of irxn --- + + // Overpotentials for the anodic and cathodic phases + scalarvalueType etaAnodic = VsV - EcorrAnodic - Phi; + + // Calculation of anodic/cathodic current density + scalarvalueType itafel = i0Anodic * exp(etaAnodic / AAnodic); + scalarvalueType iAnodic = itafel / (constV(1.0) + (itafel / iMax)); + + // Time for initial equilibration of all the interfaces present in the + // domain + + // Present simulation time + double simTime = this->currentTime; + + // Turning on the corrosion-caused motion of the interface after + // tStepsForV * timeStep time. The number 0.25 decides the width of the + // smooth transition between vMin and vMax, as described below + double arg = (simTime - tStepStartForV * userInputs.dtValue) / 0.25; + + // switchVal is zero before tStepsForV * timeStep and 1 after that time + double switchVal = 0.5 * (1.0 + tanh(arg)); + + // Normal component of the interface velocity, which switches from vMin + // to vMax at time = tStepStartForV * timeStep + + scalarvalueType vMin = constV(0.0); + scalarvalueType vMax = -constV(VMV / (zMV * FarC)) * xiAnodic * iAnodic; + scalarvalueType v = constV(1.0 - switchVal) * vMin + constV(switchVal) * vMax; + + // The mobility for all the phases + scalarvalueType MnV = MconstV * psi * abs(xiAnodic * iAnodic); + + // --- Calculation of residual terms for nAnodic --- + scalarvalueType rnAnodic = nAnodic + v * delt * magpsix; + scalargradType rnAnodicx = -MnV * delt * muAnodicx; + + // --- Calculation of residual terms for nCathodic --- + scalarvalueType rnCathodic = nCathodic; + scalargradType rnCathodicx = -MnV * delt * muCathodicx; + + // --- Calculation of residual terms for psi --- + scalarvalueType rpsi = psi - v * delt * magpsix; + scalargradType rpsix = -MnV * delt * mupsix; + + // --- Submitting the terms for the governing equations --- + // Residuals terms for the equation to evolve the order parameter + variable_list.set_scalar_value_term_RHS(0, rnAnodic); + variable_list.set_scalar_gradient_term_RHS(0, rnAnodicx); + + // Residuals terms for the equation to evolve the order parameter + variable_list.set_scalar_value_term_RHS(2, rnCathodic); + variable_list.set_scalar_gradient_term_RHS(2, rnCathodicx); + + // Residuals terms for the equation to evolve the domain parameter + variable_list.set_scalar_value_term_RHS(4, rpsi); + variable_list.set_scalar_gradient_term_RHS(4, rpsix); +} + +// ============================================================================================= +// nonExplicitEquationRHS (needed only if one or more equation is time +// independent or auxiliary) +// ============================================================================================= +// This function calculates the right-hand-side of all of the equations that are +// not explicit time-dependent equations. It takes "variable_list" as an input, +// which is a list of the value and derivatives of each of the variables at a +// specific quadrature point. The (x,y,z) location of that quadrature point is +// given by "q_point_loc". The function outputs two terms to variable_list -- +// one proportional to the test function and one proportional to the gradient of +// the test function. The index for each variable in this list corresponds to +// the index given at the top of this file. + +template +void +customPDE::nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The order parameter of the anodic phase and its gradient + scalarvalueType nAnodic = variable_list.get_scalar_value(0); + scalargradType nAnodicx = variable_list.get_scalar_gradient(0); + + // The order parameter of the cathodic phase and its gradient + scalarvalueType nCathodic = variable_list.get_scalar_value(2); + scalargradType nCathodicx = variable_list.get_scalar_gradient(2); + + // The domain parameter of the electrolyte phase and its derivatives + scalarvalueType psi = variable_list.get_scalar_value(4); + scalargradType psix = variable_list.get_scalar_gradient(4); + + // The electrolyte potential and its derivatives + scalarvalueType Phi = variable_list.get_scalar_value(6); + scalargradType Phix = variable_list.get_scalar_gradient(6); + + // The reaction current density + scalarvalueType irxn = variable_list.get_scalar_value(7); + scalarvalueType xiAnodic = variable_list.get_scalar_value(8); + + // Ensuring domain and order parameters are in the range [0 1] + psi = min(psi, constV(1.0)); + psi = max(psi, constV(0.0)); + + nAnodic = min(nAnodic, constV(1.0)); + nAnodic = max(nAnodic, constV(0.0)); + + nCathodic = min(nCathodic, constV(1.0)); + nCathodic = max(nCathodic, constV(0.0)); + + // --- Calculation of terms needed in multiple expressions --- + // Magnitude of the gradient of the domain parameter + + scalarvalueType magpsix = constV(0.0); + for (int i = 0; i < dim; i++) + { + magpsix = magpsix + psix[i] * psix[i]; + } + magpsix = sqrt(magpsix); + + // --- Calculation of residual terms for muAnodic --- + + // Derivative of bulk free energy with respect to nAnodic + scalarvalueType rmuAnodic = + (nAnodic * nAnodic * nAnodic - nAnodic + constV(2.0) * gamma * nAnodic * psi * psi + + constV(2.0) * gamma * nAnodic * nCathodic * nCathodic); + + scalargradType rmuAnodicx = epssqV * nAnodicx; + + // --- Calculation of residual terms for muCathodic --- + + // Derivative of bulk free energy with respect to nCathodic + scalarvalueType rmuCathodic = (nCathodic * nCathodic * nCathodic - nCathodic + + constV(2.0) * gamma * nCathodic * psi * psi + + constV(2.0) * gamma * nCathodic * nAnodic * nAnodic); + + scalargradType rmuCathodicx = epssqV * nCathodicx; + + // --- Calculation of residual terms for mupsi --- + + // Derivative of bulk free energy with respect to psi + scalarvalueType rmupsi = + (psi * psi * psi - psi + constV(2.0) * gamma * psi * nAnodic * nAnodic + + constV(2.0) * gamma * psi * nCathodic * nCathodic); + + scalargradType rmupsix = epssqV * psix; + + // --- Calculation of irxn --- + + // Overpotential + + scalarvalueType etaAnodic = VsV - EcorrAnodic - Phi; + scalarvalueType etaCathodic = VsV - EcorrCathodic - Phi; + + // Calculation of anodic/cathodic current density + scalarvalueType itafel = i0Anodic * exp(etaAnodic / AAnodic); + scalarvalueType iAnodic = itafel / (constV(1.0) + (itafel / iMax)); + scalarvalueType iCathodic = -i0Cathodic * exp(etaCathodic / ACathodic); + + // Calculation of the interpolation factor for the anodic phase + // The factor for the cathodic phase is calculated first, and then the + // one for the anodic phase is calculated via xiAnodic = 1 - xiCathodic + scalarvalueType xiCathodic = nCathodic / max(nAnodic + nCathodic, constV(lthresh)); + + xiAnodic = constV(1.0) - xiCathodic; + + // combination of irxn + scalarvalueType rirxn = xiAnodic * iAnodic + xiCathodic * iCathodic; + + // --- Calculation of residual terms for Phi --- + scalarvalueType rPhi = -magpsix * irxn; + scalargradType rPhix = psi * kappa * Phix; + + // --- Submitting the terms for the governing equations --- + // Residuals for the equation to calculate muAnodic + variable_list.set_scalar_value_term_RHS(1, rmuAnodic); + variable_list.set_scalar_gradient_term_RHS(1, rmuAnodicx); + + // Residuals for the equation to calculate muCathodic + variable_list.set_scalar_value_term_RHS(3, rmuCathodic); + variable_list.set_scalar_gradient_term_RHS(3, rmuCathodicx); + + // Residuals for the equation to calculate mupsi + variable_list.set_scalar_value_term_RHS(5, rmupsi); + variable_list.set_scalar_gradient_term_RHS(5, rmupsix); + + // Residuals for the equation to evolve the Potential + variable_list.set_scalar_value_term_RHS(6, rPhi); + variable_list.set_scalar_gradient_term_RHS(6, rPhix); + + // irxn for the equation to evolve irxn + variable_list.set_scalar_value_term_RHS(7, rirxn); + + // Equation to solve for xiAnodic + variable_list.set_scalar_value_term_RHS(8, xiAnodic); +} + +// ============================================================================================= +// equationLHS (needed only if at least one equation is time independent) +// ============================================================================================= +// This function calculates the left-hand-side of time-independent equations. It +// takes "variable_list" as an input, which is a list of the value and +// derivatives of each of the variables at a specific quadrature point. The +// (x,y,z) location of that quadrature point is given by "q_point_loc". The +// function outputs two terms to variable_list -- one proportional to the test +// function and one proportional to the gradient of the test function -- for the +// left-hand-side of the equation. The index for each variable in this list +// corresponds to the index given at the top of this file. If there are multiple +// elliptic equations, conditional statements should be sed to ensure that the +// correct residual is being submitted. The index of the field being solved can +// be accessed by "this->currentFieldIndex". + +template +void +customPDE::equationLHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // The order parameter of the anodic phase + scalarvalueType nAnodic = variable_list.get_scalar_value(0); + + // The order parameter of the cathodic phase + scalarvalueType nCathodic = variable_list.get_scalar_value(2); + + // The domain parameter of the electrolyt phase and its derivatives + scalarvalueType psi = variable_list.get_scalar_value(4); + scalargradType psix = variable_list.get_scalar_gradient(4); + + // The electrolyte potential and its derivatives + scalarvalueType Phi = variable_list.get_scalar_value(6); + + // The change in potential in the electrode and its derivatives + scalarvalueType DPhi = variable_list.get_change_in_scalar_value(6); + scalargradType DPhix = variable_list.get_change_in_scalar_gradient(6); + + // The interpolation factor for the anodic phase + scalarvalueType xiAnodic = variable_list.get_scalar_value(8); + + // Ensuring that the domain and order parameters remain within the range + // [0 1] + psi = min(psi, constV(1.0)); + psi = max(psi, constV(0.0)); + + nAnodic = min(nAnodic, constV(1.0)); + nAnodic = max(nAnodic, constV(0.0)); + + nCathodic = min(nCathodic, constV(1.0)); + nCathodic = max(nCathodic, constV(0.0)); + + // --- Calculation of terms needed in multiple expressions --- + + // Magnitude of the gradient of the domain parameter + scalarvalueType magpsix = constV(0.0); + for (int i = 0; i < dim; i++) + { + magpsix = magpsix + psix[i] * psix[i]; + } + magpsix = sqrt(magpsix); + + // --- Calculation of residual terms for DPhi (delta phi on the LHS) --- + // --- Calculation of residual terms for irxnp (the derivative of + // irxn wrt to Phi) --- + + // Overpotential + scalarvalueType etaAnodic = VsV - EcorrAnodic - Phi; + scalarvalueType etaCathodic = VsV - EcorrCathodic - Phi; + + // Calculation of anodic/cathodic current density + scalarvalueType itafel = i0Anodic * exp(etaAnodic / AAnodic); + scalarvalueType iCathodic = -i0Cathodic * exp(etaCathodic / ACathodic); + + // The interpolation factor for the cathodic phase + scalarvalueType xiCathodic = constV(1.0) - xiAnodic; + + // The derivative of irxn wrt to Phi + scalarvalueType irxnp = + -xiAnodic * (iMax / (iMax + itafel)) * (iMax / (iMax + itafel)) * (itafel / AAnodic) - + xiCathodic * (iCathodic / ACathodic); + + scalarvalueType rDPhi = magpsix * irxnp * DPhi; + scalargradType rDPhix = -psi * kappa * DPhix; + + // Residuals for the equation to evolve the potential(Phi) + variable_list.set_scalar_value_term_LHS(6, rDPhi); + variable_list.set_scalar_gradient_term_LHS(6, rDPhix); +} diff --git a/automatic_tests/corrosion_microgalvanic/gold_integratedFields.txt b/automatic_tests/corrosion_microgalvanic/gold_integratedFields.txt new file mode 100644 index 000000000..11ea774a7 --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/gold_integratedFields.txt @@ -0,0 +1,9 @@ +0 f_tot 49.99999772 +10 f_tot 50.27615732 +20 f_tot 50.55017018 +30 f_tot 50.81957617 +40 f_tot 51.08697849 +50 f_tot 51.35331882 +60 f_tot 51.62122273 +70 f_tot 51.89124109 +80 f_tot 52.1630525 diff --git a/automatic_tests/corrosion_microgalvanic/parameters.prm b/automatic_tests/corrosion_microgalvanic/parameters.prm new file mode 100644 index 000000000..eadf0342e --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/parameters.prm @@ -0,0 +1,165 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +# Each axes spans from zero to the specified length Unit: m +set Domain size X = 1e-6 +set Domain size Y = 1e-6 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 4 +set Subdivisions Y = 4 + +set Refine factor = 5 + +set Element degree = 1 + +# ================================================================================= +# Set the adaptive mesh refinement parameters +# ================================================================================= +set Mesh adaptivity = true + +set Max refinement level = 6 +set Min refinement level = 2 + +subsection Refinement criterion: psi + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: nAnodic + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: nCathodic + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +set Steps between remeshing operations = 2000 + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 0.01 + +set Number of time steps = 8000 + +# ================================================================================= +# Set the linear solver parameters +# ================================================================================= + +subsection Linear solver parameters: Phi + set Tolerance type = ABSOLUTE_RESIDUAL + + set Tolerance value = 1.0e-7 + + set Maximum linear solver iterations = 1000000 +end + +# ================================================================================= +# Set the nonlinear solver parameters +# ================================================================================= + +set Maximum nonlinear solver iterations = 1000 + +subsection Nonlinear solver parameters: Phi + set Tolerance type = ABSOLUTE_SOLUTION_CHANGE + set Tolerance value = 1.0e-10 + set Use backtracking line search damping = false + set Constant damping value = 0.3 + set Use Laplace's equation to determine the initial guess = true +end + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= + +set Boundary condition for variable nAnodic = NATURAL +set Boundary condition for variable muAnodic = NATURAL +set Boundary condition for variable nCathodic = NATURAL +set Boundary condition for variable muCathodic = NATURAL +set Boundary condition for variable psi = NATURAL +set Boundary condition for variable mupsi = NATURAL +set Boundary condition for variable Phi = NATURAL +set Boundary condition for variable irxn = NATURAL +set Boundary condition for variable xiAnodic = NATURAL + +# ================================================================================= +# Set the output parameters +# ================================================================================= + +set Output condition = EQUAL_SPACING + +set Number of outputs = 8 + +set Skip print steps = 1000 + +# ================================================================================= +# Set the model constants +# ================================================================================= +# For this application all units are SI + +#The dissolving metal cation molar volume (m^3/mol) +set Model constant VMV= 1.3736e-5, DOUBLE + +# The dissolving metal cation charge number +set Model constant zMV = 2.0, DOUBLE + +# The gradient energy coefficient, epsilon^2 +# (divided by the bulk energy coefficient, W) +set Model constant epssqV = 0.31e-16, DOUBLE + +# Corrosion potential for Anodic, V +set Model constant EcorrAnodic = -1.424, DOUBLE + +# Corrosion potential for Cathodic, V +set Model constant EcorrCathodic = -1.151, DOUBLE + +# Tafel slope for Anodic Unit: V +set Model constant AAnodic = 0.003474, DOUBLE + +# Tafel slope for Cathodic Unit: V +set Model constant ACathodic = -0.03626, DOUBLE + +# Corrosion current density for Anodic, A/m^2 +set Model constant i0Anodic = 0.081, DOUBLE + +# Corrosion current density for Cathodic, A/m^2 +set Model constant i0Cathodic = 0.017, DOUBLE + +# Applied potential (V) +set Model constant VsV = 0.0, DOUBLE + +# Safety factor for domain parameter to avoid division by zero +set Model constant lthresh = 1.0e-3, DOUBLE + +# Constant gamma for bulk free energy (equal to 3/2) +set Model constant gamma = 1.5, DOUBLE + +# variable KAPPA for electrolyte conductivity in comsol S/m +set Model constant kappa = 0.001, DOUBLE + +# Limitation of anodic current density A/m^2 +set Model constant iMax = 100.0, DOUBLE + +# Time step at which velocity is switched on +set Model constant tStepStartForV = 0000.0, DOUBLE + +# Size of beta particle, m +set Model constant cathodeThickness = 0.2e-6, DOUBLE + +# initial guess for Phi, V +set Model constant guessValPhi = 1.409, DOUBLE + diff --git a/automatic_tests/corrosion_microgalvanic/postprocess.cc b/automatic_tests/corrosion_microgalvanic/postprocess.cc new file mode 100644 index 000000000..928a6dae1 --- /dev/null +++ b/automatic_tests/corrosion_microgalvanic/postprocess.cc @@ -0,0 +1,57 @@ +// ================================================================================= +// Set the attributes of the postprocessing variables +// ================================================================================= +// This function is analogous to 'loadVariableAttributes' in 'equations.cc', but +// for the postprocessing expressions. It sets the attributes for each +// postprocessing expression, including its name, whether it is a vector or +// scalar (only scalars are supported at present), its dependencies on other +// variables and their derivatives, and whether to calculate an integral of the +// postprocessed quantity over the entire domain. + +void +variableAttributeLoader::loadPostProcessorVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "f_tot"); + set_variable_type(0, SCALAR); + + set_dependencies_value_term_RHS(0, "psi"); + set_dependencies_gradient_term_RHS(0, ""); + + set_output_integral(0, true); +} + +// ============================================================================================= +// postProcessedFields: Set the postprocessing expressions +// ============================================================================================= +// This function is analogous to 'explicitEquationRHS' and +// 'nonExplicitEquationRHS' in equations.cc. It takes in "variable_list" and +// "q_point_loc" as inputs and outputs two terms in the expression for the +// postprocessing variable -- one proportional to the test function and one +// proportional to the gradient of the test function. The index for each +// variable in this list corresponds to the index given at the top of this file +// (for submitting the terms) and the index in 'equations.cc' for assigning the +// values/derivatives of the primary variables. + +template +void +customPDE::postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The order parameter and its derivatives + scalarvalueType psi = variable_list.get_scalar_value(4); + + // --- Setting the expressions for the terms in the postprocessing expressions + // --- + + // --- Submitting the terms for the postprocessing expressions --- + + pp_variable_list.set_scalar_value_term_RHS(0, psi * 10e13); +} diff --git a/automatic_tests/coupledCahnHilliardAllenCahn/CMakeLists.txt b/automatic_tests/coupledCahnHilliardAllenCahn/CMakeLists.txt new file mode 100644 index 000000000..82b46ea0a --- /dev/null +++ b/automatic_tests/coupledCahnHilliardAllenCahn/CMakeLists.txt @@ -0,0 +1,165 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/ICs_and_BCs.cc b/automatic_tests/coupledCahnHilliardAllenCahn/ICs_and_BCs.cc similarity index 100% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/ICs_and_BCs.cc rename to automatic_tests/coupledCahnHilliardAllenCahn/ICs_and_BCs.cc diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/customPDE.h b/automatic_tests/coupledCahnHilliardAllenCahn/customPDE.h similarity index 100% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/customPDE.h rename to automatic_tests/coupledCahnHilliardAllenCahn/customPDE.h diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/equations.cc b/automatic_tests/coupledCahnHilliardAllenCahn/equations.cc similarity index 100% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/equations.cc rename to automatic_tests/coupledCahnHilliardAllenCahn/equations.cc diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/gold_integratedFields.txt b/automatic_tests/coupledCahnHilliardAllenCahn/gold_integratedFields.txt similarity index 100% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/gold_integratedFields.txt rename to automatic_tests/coupledCahnHilliardAllenCahn/gold_integratedFields.txt diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/parameters.prm b/automatic_tests/coupledCahnHilliardAllenCahn/parameters.prm similarity index 100% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/parameters.prm rename to automatic_tests/coupledCahnHilliardAllenCahn/parameters.prm diff --git a/tests/automatic_tests/coupledCahnHilliardAllenCahn/postprocess.cc b/automatic_tests/coupledCahnHilliardAllenCahn/postprocess.cc similarity index 100% rename from tests/automatic_tests/coupledCahnHilliardAllenCahn/postprocess.cc rename to automatic_tests/coupledCahnHilliardAllenCahn/postprocess.cc diff --git a/automatic_tests/grainGrowth/CMakeLists.txt b/automatic_tests/grainGrowth/CMakeLists.txt new file mode 100644 index 000000000..82b46ea0a --- /dev/null +++ b/automatic_tests/grainGrowth/CMakeLists.txt @@ -0,0 +1,165 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/automatic_tests/grainGrowth/ICs_and_BCs.cc b/automatic_tests/grainGrowth/ICs_and_BCs.cc new file mode 100644 index 000000000..a0db76df7 --- /dev/null +++ b/automatic_tests/grainGrowth/ICs_and_BCs.cc @@ -0,0 +1,181 @@ +// =========================================================================== +// FUNCTION FOR INITIAL CONDITIONS +// =========================================================================== + +template +void +customPDE::setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) +{ + // --------------------------------------------------------------------- + // ENTER THE INITIAL CONDITIONS HERE + // --------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the initial condition for each variable + // according to its variable index + + // The initial condition is a set of overlapping circles/spheres defined + // by a hyperbolic tangent function. The center of each circle/sphere is + // given by "center_list" and its radius is given by "radius". + + if (index < 5) + { + std::vector> center_list; + + // The big grains + { + Point center(0.2, 0.15); + center_list.push_back(center); + } + { + Point center(0.25, 0.7); + center_list.push_back(center); + } + { + Point center(0.5, 0.5); + center_list.push_back(center); + } + { + Point center(0.6, 0.85); + center_list.push_back(center); + } + { + Point center(0.85, 0.35); + center_list.push_back(center); + } + + // The medium grains + { + Point center(0.08, 0.92); + center_list.push_back(center); + } + { + Point center(0.75, 0.6); + center_list.push_back(center); + } + { + Point center(0.75, 0.1); + center_list.push_back(center); + } + { + Point center(0.2, 0.45); + center_list.push_back(center); + } + { + Point center(0.85, 0.85); + center_list.push_back(center); + } + + // The small grains + { + Point center(0.55, 0.05); + center_list.push_back(center); + } + { + Point center(0.1, 0.35); + center_list.push_back(center); + } + { + Point center(0.95, 0.65); + center_list.push_back(center); + } + { + Point center(0.9, 0.15); + center_list.push_back(center); + } + { + Point center(0.45, 0.25); + center_list.push_back(center); + } + + std::vector rad = {0.14, + 0.14, + 0.14, + 0.14, + 0.14, + 0.08, + 0.08, + 0.08, + 0.08, + 0.08, + 0.05, + 0.05, + 0.05, + 0.05, + 0.05}; + + double dist = 0.0; + scalar_IC = 0; + + for (unsigned int dir = 0; dir < dim; dir++) + { + dist += (p[dir] - center_list[index][dir] * userInputs.domain_size[dir]) * + (p[dir] - center_list[index][dir] * userInputs.domain_size[dir]); + } + dist = std::sqrt(dist); + + scalar_IC += + 0.5 * (1.0 - std::tanh((dist - rad[index] * userInputs.domain_size[0]) / 0.5)); + + dist = 0.0; + for (unsigned int dir = 0; dir < dim; dir++) + { + dist += (p[dir] - center_list[index + 5][dir] * userInputs.domain_size[dir]) * + (p[dir] - center_list[index + 5][dir] * userInputs.domain_size[dir]); + } + dist = std::sqrt(dist); + + scalar_IC += + 0.5 * + (1.0 - std::tanh((dist - rad[index + 5] * userInputs.domain_size[0]) / 0.5)); + + dist = 0.0; + for (unsigned int dir = 0; dir < dim; dir++) + { + dist += (p[dir] - center_list[index + 10][dir] * userInputs.domain_size[dir]) * + (p[dir] - center_list[index + 10][dir] * userInputs.domain_size[dir]); + } + dist = std::sqrt(dist); + + scalar_IC += + 0.5 * + (1.0 - std::tanh((dist - rad[index + 10] * userInputs.domain_size[0]) / 0.5)); + } + else + { + scalar_IC = 0.0; + } + + // -------------------------------------------------------------------------- +} + +// =========================================================================== +// FUNCTION FOR NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS +// =========================================================================== + +template +void +customPDE::setNonUniformDirichletBCs( + [[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) +{ + // -------------------------------------------------------------------------- + // ENTER THE NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS HERE + // -------------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the boundary condition for each variable + // according to its variable index. This function can be left blank if there + // are no non-uniform Dirichlet boundary conditions. For BCs that change in + // time, you can access the current time through the variable "time". The + // boundary index can be accessed via the variable "direction", which starts + // at zero and uses the same order as the BC specification in parameters.in + // (i.e. left = 0, right = 1, bottom = 2, top = 3, front = 4, back = 5). + + // ------------------------------------------------------------------------- +} diff --git a/automatic_tests/grainGrowth/customPDE.h b/automatic_tests/grainGrowth/customPDE.h new file mode 100644 index 000000000..38a1580a0 --- /dev/null +++ b/automatic_tests/grainGrowth/customPDE.h @@ -0,0 +1,94 @@ +#include + +using namespace dealii; + +template +class customPDE : public MatrixFreePDE +{ +public: + // Constructor + customPDE(userInputParameters _userInputs) + : MatrixFreePDE(_userInputs) + , userInputs(_userInputs) {}; + + // Function to set the initial conditions (in ICs_and_BCs.h) + void + setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) override; + + // Function to set the non-uniform Dirichlet boundary conditions (in + // ICs_and_BCs.h) + void + setNonUniformDirichletBCs([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) override; + +private: +#include + + const userInputParameters userInputs; + + // Function to set the RHS of the governing equations for explicit time + // dependent equations (in equations.h) + void + explicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the RHS of the governing equations for all other equations + // (in equations.h) + void + nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the LHS of the governing equations (in equations.h) + void + equationLHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + +// Function to set postprocessing expressions (in postprocess.h) +#ifdef POSTPROCESS_FILE_EXISTS + void + postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; +#endif + +// Function to set the nucleation probability (in nucleation.h) +#ifdef NUCLEATION_FILE_EXISTS + double + getNucleationProbability([[maybe_unused]] variableValueContainer variable_value, + [[maybe_unused]] double dV) const override; +#endif + + // ================================================================ + // Methods specific to this subclass + // ================================================================ + + // ================================================================ + // Model constants specific to this subclass + // ================================================================ + + double MnV = userInputs.get_model_constant_double("MnV"); + double KnV = userInputs.get_model_constant_double("KnV"); + double alpha = userInputs.get_model_constant_double("alpha"); + + // ================================================================ +}; diff --git a/automatic_tests/grainGrowth/equations.cc b/automatic_tests/grainGrowth/equations.cc new file mode 100644 index 000000000..34542af81 --- /dev/null +++ b/automatic_tests/grainGrowth/equations.cc @@ -0,0 +1,141 @@ +// ================================================================================= +// Set the attributes of the primary field variables +// ================================================================================= +// This function sets attributes for each variable/equation in the app. The +// attributes are set via standardized function calls. The first parameter for +// each function call is the variable index (starting at zero). The first set of +// variable/equation attributes are the variable name (any string), the variable +// type (SCALAR/VECTOR), and the equation type (EXPLICIT_TIME_DEPENDENT/ +// TIME_INDEPENDENT/AUXILIARY). The next set of attributes describe the +// dependencies for the governing equation on the values and derivatives of the +// other variables for the value term and gradient term of the RHS and the LHS. +// The final pair of attributes determine whether a variable represents a field +// that can nucleate and whether the value of the field is needed for nucleation +// rate calculations. + +void +variableAttributeLoader::loadVariableAttributes() +{ + for (unsigned int var_index = 0; var_index < 6; var_index++) + { + std::string var_name = "n"; + var_name.append(std::to_string(var_index)); + + set_variable_name(var_index, var_name); + set_variable_type(var_index, SCALAR); + set_variable_equation_type(var_index, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(var_index, "n0, n1, n2, n3, n4, n5"); + set_dependencies_gradient_term_RHS( + var_index, + "grad(n0), grad(n1), grad(n2), grad(n3), grad(n4), grad(n5)"); + } +} + +// ============================================================================================= +// explicitEquationRHS (needed only if one or more equation is explict time +// dependent) +// ============================================================================================= +// This function calculates the right-hand-side of the explicit time-dependent +// equations for each variable. It takes "variable_list" as an input, which is a +// list of the value and derivatives of each of the variables at a specific +// quadrature point. The (x,y,z) location of that quadrature point is given by +// "q_point_loc". The function outputs two terms to variable_list -- one +// proportional to the test function and one proportional to the gradient of the +// test function. The index for each variable in this list corresponds to the +// index given at the top of this file. + +template +void +customPDE::explicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + VectorizedArray fnV = constV(0.0); + scalarvalueType ni, nj; + scalargradType nix; + + // In this application, create temporary variables for the residual terms. We + // cannot call 'set_scalar_value_residual_term' and + // 'set_scalar_gradient_residual_term'in the for loop below because those + // functions write over the scalar value and scalar gradient internal + // variables in 'variable_list' (for performance reasons). Therefore, we wait + // to set the residual terms until all the residuals have been calculated. + + std::vector value_terms; + value_terms.resize(userInputs.number_of_variables); + std::vector gradient_terms; + gradient_terms.resize(userInputs.number_of_variables); + + for (unsigned int i = 0; i < userInputs.number_of_variables; i++) + { + ni = variable_list.get_scalar_value(i); + nix = variable_list.get_scalar_gradient(i); + fnV = -ni + ni * ni * ni; + for (unsigned int j = 0; j < userInputs.number_of_variables; j++) + { + if (i != j) + { + nj = variable_list.get_scalar_value(j); + fnV += constV(2.0 * alpha) * ni * nj * nj; + } + } + value_terms[i] = ni - constV(userInputs.dtValue * MnV) * fnV; + gradient_terms[i] = constV(-userInputs.dtValue * KnV * MnV) * nix; + } + + // --- Submitting the terms for the governing equations --- + + for (unsigned int i = 0; i < userInputs.number_of_variables; i++) + { + variable_list.set_scalar_value_term_RHS(i, value_terms[i]); + variable_list.set_scalar_gradient_term_RHS(i, gradient_terms[i]); + } +} + +// ============================================================================================= +// nonExplicitEquationRHS (needed only if one or more equation is time +// independent or auxiliary) +// ============================================================================================= +// This function calculates the right-hand-side of all of the equations that are +// not explicit time-dependent equations. It takes "variable_list" as an input, +// which is a list of the value and derivatives of each of the variables at a +// specific quadrature point. The (x,y,z) location of that quadrature point is +// given by "q_point_loc". The function outputs two terms to variable_list -- +// one proportional to the test function and one proportional to the gradient of +// the test function. The index for each variable in this list corresponds to +// the index given at the top of this file. + +template +void +customPDE::nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{} + +// ============================================================================================= +// equationLHS (needed only if at least one equation is time independent) +// ============================================================================================= +// This function calculates the left-hand-side of time-independent equations. It +// takes "variable_list" as an input, which is a list of the value and +// derivatives of each of the variables at a specific quadrature point. The +// (x,y,z) location of that quadrature point is given by "q_point_loc". The +// function outputs two terms to variable_list -- one proportional to the test +// function and one proportional to the gradient of the test function -- for the +// left-hand-side of the equation. The index for each variable in this list +// corresponds to the index given at the top of this file. If there are multiple +// elliptic equations, conditional statements should be sed to ensure that the +// correct residual is being submitted. The index of the field being solved can +// be accessed by "this->currentFieldIndex". + +template +void +customPDE::equationLHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{} diff --git a/automatic_tests/grainGrowth/gold_integratedFields.txt b/automatic_tests/grainGrowth/gold_integratedFields.txt new file mode 100644 index 000000000..7f60e4cba --- /dev/null +++ b/automatic_tests/grainGrowth/gold_integratedFields.txt @@ -0,0 +1,21 @@ +0 f_tot 6487.3288 +2 f_tot 9523.096909 +4 f_tot 12228.01679 +6 f_tot 14571.46539 +8 f_tot 16353.94815 +10 f_tot 27053.87822 +12 f_tot 28246.73123 +14 f_tot 29073.71992 +16 f_tot 29639.01096 +18 f_tot 30060.30424 +20 f_tot 29746.59861 +22 f_tot 29906.94022 +24 f_tot 30009.04225 +26 f_tot 30064.36572 +28 f_tot 30081.18429 +30 f_tot 30086.33837 +32 f_tot 30089.74429 +34 f_tot 30096.49583 +36 f_tot 30101.63484 +38 f_tot 30107.48216 +40 f_tot 30113.19384 diff --git a/automatic_tests/grainGrowth/parameters.prm b/automatic_tests/grainGrowth/parameters.prm new file mode 100644 index 000000000..19ddf14e7 --- /dev/null +++ b/automatic_tests/grainGrowth/parameters.prm @@ -0,0 +1,119 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +set Domain size X = 100 +set Domain size Y = 100 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 3 +set Subdivisions Y = 3 +set Subdivisions Z = 3 + +set Refine factor = 4 + +set Element degree = 2 + +# ================================================================================= +# Set the adaptive mesh refinement parameters +# ================================================================================= +set Mesh adaptivity = true + +set Max refinement level = 6 +set Min refinement level = 4 + +set Steps between remeshing operations = 20 + +subsection Refinement criterion: n0 + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: n1 + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: n2 + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: n3 + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: n4 + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +subsection Refinement criterion: n5 + set Criterion type = VALUE + set Value lower bound = 1e-3 + set Value upper bound = 0.999 +end + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 1.0e-1 + +set Number of time steps = 400 + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 20 + +set Skip print steps = 50 + +# ================================================================================= +# Set the grain reassignment parameters +# ================================================================================= + +set Activate grain reassignment = true +set Order parameter fields for grain reassignment = n0, n1, n2, n3, n4, n5 +set Time steps between grain reassignments = 100 +set Order parameter cutoff for grain identification = 1e-2 +set Buffer between grains before reassignment = 5.0 +set Load grain structure = false + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable n0 = NATURAL +set Boundary condition for variable n1 = NATURAL +set Boundary condition for variable n2 = NATURAL +set Boundary condition for variable n3 = NATURAL +set Boundary condition for variable n4 = NATURAL +set Boundary condition for variable n5 = NATURAL + +# ================================================================================= +# Set the model constants +# ================================================================================= + +# The AC mobility, MnV in equations.h +set Model constant MnV = 1.0, DOUBLE + +# Gradient energy coefficient +set Model constant KnV = 0.1, DOUBLE + +# Grain interaction coeffiecient +set Model constant alpha = 1.5, DOUBLE diff --git a/automatic_tests/grainGrowth/postprocess.cc b/automatic_tests/grainGrowth/postprocess.cc new file mode 100644 index 000000000..f64d41994 --- /dev/null +++ b/automatic_tests/grainGrowth/postprocess.cc @@ -0,0 +1,130 @@ +// ================================================================================= +// Set the attributes of the postprocessing variables +// ================================================================================= +// This function is analogous to 'loadVariableAttributes' in 'equations.h', but +// for the postprocessing expressions. It sets the attributes for each +// postprocessing expression, including its name, whether it is a vector or +// scalar (only scalars are supported at present), its dependencies on other +// variables and their derivatives, and whether to calculate an integral of the +// postprocessed quantity over the entire domain. + +void +variableAttributeLoader::loadPostProcessorVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "feature_ids"); + set_variable_type(0, SCALAR); + + set_dependencies_value_term_RHS(0, "n0, n1, n2, n3, n4, n5"); + set_dependencies_gradient_term_RHS(0, ""); + + set_output_integral(0, false); + + // Variable 1 + set_variable_name(1, "f_tot"); + set_variable_type(1, SCALAR); + + set_dependencies_value_term_RHS(1, "n0, n1, n2, n3, n4, n5"); + set_dependencies_gradient_term_RHS(1, ""); + + set_output_integral(1, true); +} + +// ============================================================================================= +// postProcessedFields: Set the postprocessing expressions +// ============================================================================================= +// This function is analogous to 'explicitEquationRHS' and +// 'nonExplicitEquationRHS' in equations.h. It takes in "variable_list" and +// "q_point_loc" as inputs and outputs two terms in the expression for the +// postprocessing variable -- one proportional to the test function and one +// proportional to the gradient of the test function. The index for each +// variable in this list corresponds to the index given at the top of this file +// (for submitting the terms) and the index in 'equations.h' for assigning the +// values/derivatives of the primary variables. + +template +void +customPDE::postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + scalarvalueType ni; + + scalarvalueType max_val = constV(-1.0); + scalarvalueType max_op = constV(100.0); + for (unsigned int i = 0; i < userInputs.number_of_variables; i++) + { + ni = variable_list.get_scalar_value(i); + + for (unsigned int v = 0; v < ni.size(); v++) + { + if (ni[v] > max_val[v]) + { + max_val[v] = ni[v]; + max_op[v] = i; + } + } + } + + scalarvalueType feature_ids = constV(-1.0); + for (unsigned int v = 0; v < ni.size(); v++) + { + for (unsigned int g = 0; g < this->simplified_grain_representations.size(); g++) + { + unsigned int max_op_nonvec = (unsigned int) std::abs(max_op[v]); + + if (this->simplified_grain_representations[g].getOrderParameterId() == + max_op_nonvec) + { + Point q_point_loc_nonvec; + for (unsigned int d = 0; d < dim; d++) + { + q_point_loc_nonvec(d) = q_point_loc(d)[v]; + } + + double dist = 0.0; + for (unsigned int d = 0; d < dim; d++) + { + dist += (q_point_loc_nonvec(d) - + this->simplified_grain_representations[g].getCenter()(d)) * + (q_point_loc_nonvec(d) - + this->simplified_grain_representations[g].getCenter()(d)); + } + dist = std::sqrt(dist); + + if (dist < (this->simplified_grain_representations[g].getRadius() + + userInputs.buffer_between_grains / 2.0)) + { + feature_ids[v] = + (double) (this->simplified_grain_representations[g].getGrainId()); + } + } + } + } + + scalarvalueType sum_n = constV(0.0); + for (unsigned int i = 0; i < userInputs.number_of_variables; i++) + { + ni = variable_list.get_scalar_value(i); + sum_n += ni; + } + for (unsigned int v = 0; v < ni.size(); v++) + { + if (sum_n[v] < 0.01) + { + max_op[v] = -1.0; + feature_ids[v] = -1.0; + } + } + + // --- Submitting the terms for the postprocessing expressions --- + + pp_variable_list.set_scalar_value_term_RHS(0, feature_ids); + pp_variable_list.set_scalar_value_term_RHS(1, max_op); +} diff --git a/tests/automatic_tests/main.cc b/automatic_tests/main.cc similarity index 100% rename from tests/automatic_tests/main.cc rename to automatic_tests/main.cc diff --git a/automatic_tests/precipitateEvolution/CMakeLists.txt b/automatic_tests/precipitateEvolution/CMakeLists.txt new file mode 100644 index 000000000..82b46ea0a --- /dev/null +++ b/automatic_tests/precipitateEvolution/CMakeLists.txt @@ -0,0 +1,165 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/tests/automatic_tests/precipitateEvolution/ICs_and_BCs.cc b/automatic_tests/precipitateEvolution/ICs_and_BCs.cc similarity index 100% rename from tests/automatic_tests/precipitateEvolution/ICs_and_BCs.cc rename to automatic_tests/precipitateEvolution/ICs_and_BCs.cc diff --git a/tests/automatic_tests/precipitateEvolution/customPDE.h b/automatic_tests/precipitateEvolution/customPDE.h similarity index 100% rename from tests/automatic_tests/precipitateEvolution/customPDE.h rename to automatic_tests/precipitateEvolution/customPDE.h diff --git a/tests/automatic_tests/precipitateEvolution/equations.cc b/automatic_tests/precipitateEvolution/equations.cc similarity index 100% rename from tests/automatic_tests/precipitateEvolution/equations.cc rename to automatic_tests/precipitateEvolution/equations.cc diff --git a/tests/automatic_tests/precipitateEvolution/gold_integratedFields.txt b/automatic_tests/precipitateEvolution/gold_integratedFields.txt similarity index 100% rename from tests/automatic_tests/precipitateEvolution/gold_integratedFields.txt rename to automatic_tests/precipitateEvolution/gold_integratedFields.txt diff --git a/tests/automatic_tests/precipitateEvolution/parameters.prm b/automatic_tests/precipitateEvolution/parameters.prm similarity index 100% rename from tests/automatic_tests/precipitateEvolution/parameters.prm rename to automatic_tests/precipitateEvolution/parameters.prm diff --git a/tests/automatic_tests/precipitateEvolution/postprocess.cc b/automatic_tests/precipitateEvolution/postprocess.cc similarity index 100% rename from tests/automatic_tests/precipitateEvolution/postprocess.cc rename to automatic_tests/precipitateEvolution/postprocess.cc diff --git a/automatic_tests/precipitateEvolution_pfunction/CMakeLists.txt b/automatic_tests/precipitateEvolution_pfunction/CMakeLists.txt new file mode 100644 index 000000000..82b46ea0a --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/CMakeLists.txt @@ -0,0 +1,165 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/automatic_tests/precipitateEvolution_pfunction/ICs_and_BCs.cc b/automatic_tests/precipitateEvolution_pfunction/ICs_and_BCs.cc new file mode 100644 index 000000000..bfb3b7a84 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/ICs_and_BCs.cc @@ -0,0 +1,94 @@ +// =========================================================================== +// FUNCTION FOR INITIAL CONDITIONS +// =========================================================================== + +template +void +customPDE::setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) +{ + // --------------------------------------------------------------------- + // ENTER THE INITIAL CONDITIONS HERE + // --------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the initial condition for each variable + // according to its variable index + + double center[4][3] = { + {1.0 / 3.0, 1.0 / 3.0, 0.5}, + {2.0 / 3.0, 2.0 / 3.0, 0.5}, + {3.0 / 4.0, 1.0 / 4.0, 0.5}, + {1.0 / 4.0, 3.0 / 4, 0.5} + }; + double rad[4] = {userInputs.domain_size[0] / 16.0, + userInputs.domain_size[0] / 16.0, + userInputs.domain_size[0] / 16.0, + userInputs.domain_size[0] / 16.0}; + double orientation[4] = {1, 1, 2, 3}; + double dx = userInputs.domain_size[0] / ((double) userInputs.subdivisions[0]) / + std::pow(2.0, userInputs.refine_factor); + double dist; + scalar_IC = 0; + + if (index == 0) + { + scalar_IC = 0.04; + } + + for (unsigned int i = 0; i < 4; i++) + { + dist = 0.0; + for (unsigned int dir = 0; dir < dim; dir++) + { + dist += (p[dir] - center[i][dir] * userInputs.domain_size[dir]) * + (p[dir] - center[i][dir] * userInputs.domain_size[dir]); + } + dist = std::sqrt(dist); + + if (index == orientation[i]) + { + scalar_IC += 0.5 * (1.0 - std::tanh((dist - rad[i]) / (dx))); + } + } + + if (index == 4) + { + for (unsigned int d = 0; d < dim; d++) + { + vector_IC(d) = 0.0; + } + } + + // -------------------------------------------------------------------------- +} + +// =========================================================================== +// FUNCTION FOR NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS +// =========================================================================== + +template +void +customPDE::setNonUniformDirichletBCs( + [[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) +{ + // -------------------------------------------------------------------------- + // ENTER THE NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS HERE + // -------------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the boundary condition for each variable + // according to its variable index. This function can be left blank if there + // are no non-uniform Dirichlet boundary conditions. For BCs that change in + // time, you can access the current time through the variable "time". The + // boundary index can be accessed via the variable "direction", which starts + // at zero and uses the same order as the BC specification in parameters.in + // (i.e. left = 0, right = 1, bottom = 2, top = 3, front = 4, back = 5). + + // ------------------------------------------------------------------------- +} diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/.PLibrary.hh.swp b/automatic_tests/precipitateEvolution_pfunction/PLibrary/.PLibrary.hh.swp new file mode 100644 index 000000000..eee036601 Binary files /dev/null and b/automatic_tests/precipitateEvolution_pfunction/PLibrary/.PLibrary.hh.swp differ diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/Mc.txt b/automatic_tests/precipitateEvolution_pfunction/PLibrary/Mc.txt new file mode 100644 index 000000000..d3827e75a --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/Mc.txt @@ -0,0 +1 @@ +1.0 diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/Mn.txt b/automatic_tests/precipitateEvolution_pfunction/PLibrary/Mn.txt new file mode 100644 index 000000000..bd9e0289f --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/Mn.txt @@ -0,0 +1,3 @@ +100.0 +100.0 +100.0 diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/PLibrary.cc b/automatic_tests/precipitateEvolution_pfunction/PLibrary/PLibrary.cc new file mode 100644 index 000000000..5f03c1aff --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/PLibrary.cc @@ -0,0 +1,217 @@ +// created: 2017-1-11 12:28:14 +// version: issue-1 +// url: https://github.com/prisms-center/IntegrationTools +// commit: 34f71553d3b1b3cdbbc6a02cfb8be793a16396b5 + +#ifndef PLIBRARY_CC +#define PLIBRARY_CC + +#include "PLibrary.hh" + +#include "pfunct_McV.hh" +#include "pfunct_Mn1V.hh" +#include "pfunct_Mn2V.hh" +#include "pfunct_Mn3V.hh" +#include "pfunct_faV.hh" +#include "pfunct_fbV.hh" + +#include +#include + +namespace PRISMS +{ + + void + PLibrary::checkout(std::string name, PSimpleFunction &simplefunc) + { + typedef PSimpleFunction psf; + if (name == "pfunct_faV_f") + { + simplefunc = psf(pfunct_faV_f()); + return; + } + if (name == "pfunct_faV_grad_0") + { + simplefunc = psf(pfunct_faV_grad_0()); + return; + } + if (name == "pfunct_faV_hess_0_0") + { + simplefunc = psf(pfunct_faV_hess_0_0()); + return; + } + if (name == "pfunct_fbV_f") + { + simplefunc = psf(pfunct_fbV_f()); + return; + } + if (name == "pfunct_fbV_grad_0") + { + simplefunc = psf(pfunct_fbV_grad_0()); + return; + } + if (name == "pfunct_fbV_hess_0_0") + { + simplefunc = psf(pfunct_fbV_hess_0_0()); + return; + } + if (name == "pfunct_McV_f") + { + simplefunc = psf(pfunct_McV_f()); + return; + } + if (name == "pfunct_Mn1V_f") + { + simplefunc = psf(pfunct_Mn1V_f()); + return; + } + if (name == "pfunct_Mn2V_f") + { + simplefunc = psf(pfunct_Mn2V_f()); + return; + } + if (name == "pfunct_Mn3V_f") + { + simplefunc = psf(pfunct_Mn3V_f()); + return; + } + else + throw std::runtime_error("PSimpleFunction< double*, double > " + name + + " was not found in the PLibrary"); + } + + void + PLibrary::checkout(std::string name, PFunction &func) + { + typedef PFunction pf; + if (name == "pfunct_faV") + { + func = pf(pfunct_faV()); + return; + } + if (name == "pfunct_fbV") + { + func = pf(pfunct_fbV()); + return; + } + if (name == "pfunct_McV") + { + func = pf(pfunct_McV()); + return; + } + if (name == "pfunct_Mn1V") + { + func = pf(pfunct_Mn1V()); + return; + } + if (name == "pfunct_Mn2V") + { + func = pf(pfunct_Mn2V()); + return; + } + if (name == "pfunct_Mn3V") + { + func = pf(pfunct_Mn3V()); + return; + } + throw std::runtime_error("PFunction< double*, double > " + name + + " was not found in the PLibrary"); + } + + void + PLibrary::checkout(std::string name, PSimpleBase *&simplefunc) + { + if (name == "pfunct_faV_f") + { + simplefunc = new pfunct_faV_f(); + return; + } + if (name == "pfunct_faV_grad_0") + { + simplefunc = new pfunct_faV_grad_0(); + return; + } + if (name == "pfunct_faV_hess_0_0") + { + simplefunc = new pfunct_faV_hess_0_0(); + return; + } + if (name == "pfunct_fbV_f") + { + simplefunc = new pfunct_fbV_f(); + return; + } + if (name == "pfunct_fbV_grad_0") + { + simplefunc = new pfunct_fbV_grad_0(); + return; + } + if (name == "pfunct_fbV_hess_0_0") + { + simplefunc = new pfunct_fbV_hess_0_0(); + return; + } + if (name == "pfunct_McV_f") + { + simplefunc = new pfunct_McV_f(); + return; + } + if (name == "pfunct_Mn1V_f") + { + simplefunc = new pfunct_Mn1V_f(); + return; + } + if (name == "pfunct_Mn2V_f") + { + simplefunc = new pfunct_Mn2V_f(); + return; + } + if (name == "pfunct_Mn3V_f") + { + simplefunc = new pfunct_Mn3V_f(); + return; + } + throw std::runtime_error("PSimpleBase< double*, double > " + name + + " was not found in the PLibrary"); + } + + void + PLibrary::checkout(std::string name, PFuncBase *&func) + { + if (name == "pfunct_faV") + { + func = new pfunct_faV(); + return; + } + if (name == "pfunct_fbV") + { + func = new pfunct_fbV(); + return; + } + if (name == "pfunct_McV") + { + func = new pfunct_McV(); + return; + } + if (name == "pfunct_Mn1V") + { + func = new pfunct_Mn1V(); + return; + } + if (name == "pfunct_Mn2V") + { + func = new pfunct_Mn2V(); + return; + } + if (name == "pfunct_Mn3V") + { + func = new pfunct_Mn3V(); + return; + } + throw std::runtime_error("PFuncBase< double*, double > " + name + + " was not found in the PLibrary"); + } + +} // namespace PRISMS + +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/PLibrary.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/PLibrary.hh new file mode 100644 index 000000000..488e936cc --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/PLibrary.hh @@ -0,0 +1,40 @@ +// created: 2017-1-11 12:28:14 +// version: issue-1 +// url: https://github.com/prisms-center/IntegrationTools +// commit: 34f71553d3b1b3cdbbc6a02cfb8be793a16396b5 + +#ifndef PLIBRARY_HH +#define PLIBRARY_HH + +#include +#include +#include + +namespace PRISMS +{ + + /// Library where you can find functions and basis sets + /// + namespace PLibrary + { + // Use these functions to checkout objects which manage their own memory + + void + checkout(std::string name, PSimpleFunction &simplefunc); + + void + checkout(std::string name, PFunction &func); + + // Use these functions to checkout new 'Base' objects which the user must delete + + void + checkout(std::string name, PSimpleBase *&simplefunc); + + void + checkout(std::string name, PFuncBase *&func); + + } // namespace PLibrary + +} // namespace PRISMS + +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/createPLib.py b/automatic_tests/precipitateEvolution_pfunction/PLibrary/createPLib.py new file mode 100644 index 000000000..2c668971d --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/createPLib.py @@ -0,0 +1,90 @@ +import os.path +import subprocess +import shutil +import glob +#import math +import os.path +import sys + +# ----------------------------------------------------------------------------------------- +# Function that creates a PFunction using the IntegrationToolsWriter command line interface +# ----------------------------------------------------------------------------------------- +# Inputs: +# f_name = The name the function is referred to in the code that calls the PLibrary +# var = The variable(s) used in the function express, should match the variable names in the code that calls the PLibrary +# expression = The expression for the desired function in terms of "var" +# need_grad = Flag set to true if the gradient of the expression is needed, false if not +# need_hess = Flag set to true if the hessian of the expression is needed, false if not +# dir = Directory to write the PFunction to + +def write_pfunction(f_name, var, var_name, expression, need_grad, need_hess, dir): + f_writer_string = 'fw -n ' + f_name + ' -v ' + var + ' -d ' + var_name + ' --sym "'+ expression +'"' + if need_grad == True: + f_writer_string += " --grad " + if need_hess == True: + f_writer_string += " --hess " + + f_writer_string += ' -l '+dir + + subprocess.call([f_writer_string],shell=True) + +# ----------------------------------------------------------------------------------------- +# Function that creates a PLibrary using the IntegrationTools command line interface +# ----------------------------------------------------------------------------------------- +# Inputs: +# data_type = Data type for the container holding the PLibrary +# pfunction_dir = Directory containing the PFunctions to be assembled into the PLibrary +# plibrary_dir = Directory where the PLibrary should be written + +def write_plibrary(data_type, pfunction_dir, plibrary_dir): + #l_writer_string = 'lw -d ' + pfunction_dir + ' -v "'+ data_type + '" -l '+ plibrary_dir +' -c --include ""' + l_writer_string = 'lw -d ' + pfunction_dir + ' -l ' + plibrary_dir + + + + print l_writer_string + + subprocess.call([l_writer_string],shell=True) + +# ----------------------------------------------------------------------------------------- +# Main script to generate a PLibrary +# ----------------------------------------------------------------------------------------- + +# Get parameters that will be placed in the library (eventually through the MC API, currently in text files) +coeff_file = open('fa.txt','r') +fa_coeffs = coeff_file.read().splitlines() +coeff_file.close() + +coeff_file = open('fb.txt','r') +fb_coeffs = coeff_file.read().splitlines() +coeff_file.close() + +CH_mobility_file = open('Mc.txt','r') +Mc = CH_mobility_file.read() +CH_mobility_file.close() + +AC_mobility_file = open('Mn.txt','r') +Mn = AC_mobility_file.read().splitlines() +AC_mobility_file.close() + +# Get the current directory, which is where the PFunctions and PLibrary will be created +dir = os.getcwd() + +# Write PFunctions for the free energies and their first two derivatives +write_pfunction("pfunct_faV", "c", "concentration", fa_coeffs[0]+'*c^4 +'+fa_coeffs[1]+'*c^3 + '+fa_coeffs[2]+'*c^2 +'+fa_coeffs[3]+'*c +'+fa_coeffs[4], True, True, dir) +write_pfunction("pfunct_fbV", "c", "concentration", fb_coeffs[0]+'*c^2 +'+fb_coeffs[1]+'*c +'+fb_coeffs[2], True, True, dir) + +# Write a PFunction for the Cahn-Hilliard mobility +write_pfunction("pfunct_McV", "c", "concentration", Mc, False, False, dir) + +# Write PFunctions for the Allen-Cahn mobilities +write_pfunction("pfunct_Mn1V", "n1", "concentration", Mn[0], False, False, dir) +write_pfunction("pfunct_Mn2V", "n2", "concentration", Mn[1], False, False, dir) +write_pfunction("pfunct_Mn3V", "n3", "concentration", Mn[2], False, False, dir) + +# Write the PLibrary +#write_plibrary("dealii::VectorizedArray", dir, dir) +write_plibrary("double", dir, dir) + + + diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/fa.txt b/automatic_tests/precipitateEvolution_pfunction/PLibrary/fa.txt new file mode 100644 index 000000000..9630fe8d3 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/fa.txt @@ -0,0 +1,5 @@ +1.3687 +-2.7375 +5.1622 +-4.776 +-1.6704 diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/fb.txt b/automatic_tests/precipitateEvolution_pfunction/PLibrary/fb.txt new file mode 100644 index 000000000..5617ed23c --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/fb.txt @@ -0,0 +1,3 @@ +5.0 +-5.9746 +-1.5924 diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh new file mode 100644 index 000000000..24892f24c --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_McV.hh @@ -0,0 +1,141 @@ +// created: 2017-1-11 12:28:14 +// version: master +// url: https://github.com/bpuchala/IntegrationToolsWriter.git +// commit: 13e063c3ac8e8911a726a243fdbd68f291cc58cc + +#ifndef pfunct_McV_HH +#define pfunct_McV_HH + +#include +#include +#include + +namespace PRISMS +{ + template + class pfunct_McV_f : public PSimpleBase + { + double + eval([[maybe_unused]] const VarContainer &var) const override + { + return 1.0000000000000000e+00; + } + + public: + pfunct_McV_f() + { + this->_name = "pfunct_McV_f"; + } + + std::string + csrc() const override + { + return "1.0000000000000000e+00"; + } + + std::string + sym() const override + { + return "1.0"; + } + + std::string + latex() const override + { + return "1.0"; + } + + pfunct_McV_f * + clone() const override + { + return new pfunct_McV_f(*this); + } + }; + + template + class pfunct_McV : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_McV() + { + construct(); + } + + pfunct_McV(const pfunct_McV &RHS) + : PFuncBase(RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_McV & + operator=(pfunct_McV RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_McV() + { + delete _val; + } + + pfunct_McV * + clone() const override + { + return new pfunct_McV(*this); + } + + PSimpleFunction + simplefunction() const override + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) override + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) override + { + (*_val)(var); + } + + double + operator()() const override + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_McV"; + this->_var_name.clear(); + this->_var_name.push_back("c"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_McV_f(); + } + }; + +} // namespace PRISMS +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh new file mode 100644 index 000000000..05e6dd063 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn1V.hh @@ -0,0 +1,141 @@ +// created: 2017-1-11 12:28:14 +// version: master +// url: https://github.com/bpuchala/IntegrationToolsWriter.git +// commit: 13e063c3ac8e8911a726a243fdbd68f291cc58cc + +#ifndef pfunct_Mn1V_HH +#define pfunct_Mn1V_HH + +#include +#include +#include + +namespace PRISMS +{ + template + class pfunct_Mn1V_f : public PSimpleBase + { + double + eval([[maybe_unused]] const VarContainer &var) const override + { + return 1.0000000000000000e+02; + } + + public: + pfunct_Mn1V_f() + { + this->_name = "pfunct_Mn1V_f"; + } + + std::string + csrc() const override + { + return "1.0000000000000000e+02"; + } + + std::string + sym() const override + { + return "100.0"; + } + + std::string + latex() const override + { + return "100.0"; + } + + pfunct_Mn1V_f * + clone() const override + { + return new pfunct_Mn1V_f(*this); + } + }; + + template + class pfunct_Mn1V : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_Mn1V() + { + construct(); + } + + pfunct_Mn1V(const pfunct_Mn1V &RHS) + : PFuncBase(RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_Mn1V & + operator=(pfunct_Mn1V RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_Mn1V() + { + delete _val; + } + + pfunct_Mn1V * + clone() const override + { + return new pfunct_Mn1V(*this); + } + + PSimpleFunction + simplefunction() const override + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) override + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) override + { + (*_val)(var); + } + + double + operator()() const override + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_Mn1V"; + this->_var_name.clear(); + this->_var_name.push_back("n1"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_Mn1V_f(); + } + }; + +} // namespace PRISMS +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh new file mode 100644 index 000000000..7439b6093 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn2V.hh @@ -0,0 +1,141 @@ +// created: 2017-1-11 12:28:14 +// version: master +// url: https://github.com/bpuchala/IntegrationToolsWriter.git +// commit: 13e063c3ac8e8911a726a243fdbd68f291cc58cc + +#ifndef pfunct_Mn2V_HH +#define pfunct_Mn2V_HH + +#include +#include +#include + +namespace PRISMS +{ + template + class pfunct_Mn2V_f : public PSimpleBase + { + double + eval([[maybe_unused]] const VarContainer &var) const override + { + return 1.0000000000000000e+02; + } + + public: + pfunct_Mn2V_f() + { + this->_name = "pfunct_Mn2V_f"; + } + + std::string + csrc() const override + { + return "1.0000000000000000e+02"; + } + + std::string + sym() const override + { + return "100.0"; + } + + std::string + latex() const override + { + return "100.0"; + } + + pfunct_Mn2V_f * + clone() const override + { + return new pfunct_Mn2V_f(*this); + } + }; + + template + class pfunct_Mn2V : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_Mn2V() + { + construct(); + } + + pfunct_Mn2V(const pfunct_Mn2V &RHS) + : PFuncBase(RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_Mn2V & + operator=(pfunct_Mn2V RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_Mn2V() + { + delete _val; + } + + pfunct_Mn2V * + clone() const override + { + return new pfunct_Mn2V(*this); + } + + PSimpleFunction + simplefunction() const override + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) override + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) override + { + (*_val)(var); + } + + double + operator()() const override + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_Mn2V"; + this->_var_name.clear(); + this->_var_name.push_back("n2"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_Mn2V_f(); + } + }; + +} // namespace PRISMS +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh new file mode 100644 index 000000000..ae765e522 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_Mn3V.hh @@ -0,0 +1,141 @@ +// created: 2017-1-11 12:28:14 +// version: master +// url: https://github.com/bpuchala/IntegrationToolsWriter.git +// commit: 13e063c3ac8e8911a726a243fdbd68f291cc58cc + +#ifndef pfunct_Mn3V_HH +#define pfunct_Mn3V_HH + +#include +#include +#include + +namespace PRISMS +{ + template + class pfunct_Mn3V_f : public PSimpleBase + { + double + eval([[maybe_unused]] const VarContainer &var) const override + { + return 1.0000000000000000e+02; + } + + public: + pfunct_Mn3V_f() + { + this->_name = "pfunct_Mn3V_f"; + } + + std::string + csrc() const override + { + return "1.0000000000000000e+02"; + } + + std::string + sym() const override + { + return "100.0"; + } + + std::string + latex() const override + { + return "100.0"; + } + + pfunct_Mn3V_f * + clone() const override + { + return new pfunct_Mn3V_f(*this); + } + }; + + template + class pfunct_Mn3V : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_Mn3V() + { + construct(); + } + + pfunct_Mn3V(const pfunct_Mn3V &RHS) + : PFuncBase(RHS) + { + construct(false); + + _val = RHS._val->clone(); + } + + pfunct_Mn3V & + operator=(pfunct_Mn3V RHS) + { + using std::swap; + + swap(_val, RHS._val); + + return *this; + } + + ~pfunct_Mn3V() + { + delete _val; + } + + pfunct_Mn3V * + clone() const override + { + return new pfunct_Mn3V(*this); + } + + PSimpleFunction + simplefunction() const override + { + return PSimpleFunction(*_val); + } + + double + operator()(const VarContainer &var) override + { + return (*_val)(var); + } + + void + eval(const VarContainer &var) override + { + (*_val)(var); + } + + double + operator()() const override + { + return (*_val)(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_Mn3V"; + this->_var_name.clear(); + this->_var_name.push_back("n3"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + if (!allocate) + return; + + _val = new pfunct_Mn3V_f(); + } + }; + +} // namespace PRISMS +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh new file mode 100644 index 000000000..8b9d14008 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_faV.hh @@ -0,0 +1,300 @@ +// created: 2017-1-11 12:28:13 +// version: master +// url: https://github.com/bpuchala/IntegrationToolsWriter.git +// commit: 13e063c3ac8e8911a726a243fdbd68f291cc58cc + +#ifndef pfunct_faV_HH +#define pfunct_faV_HH + +#include +#include +#include + +namespace PRISMS +{ + template + class pfunct_faV_f : public PSimpleBase + { + double + eval(const VarContainer &var) const override + { + return 5.1622000000000003e+00 * (var[0] * var[0]) + + -2.7374999999999998e+00 * (var[0] * var[0] * var[0]) + + -4.7759999999999998e+00 * var[0] + + 1.3687000000000000e+00 * ((var[0] * var[0]) * (var[0] * var[0])) - + 1.6704000000000001e+00; + } + + public: + pfunct_faV_f() + { + this->_name = "pfunct_faV_f"; + } + + std::string + csrc() const override + { + return " 5.1622000000000003e+00*(var[0]*var[0])+-2.7374999999999998e+00*(var[0]*" + "var[0]*var[0])+-4.7759999999999998e+00*var[0]+1.3687000000000000e+00*((var[" + "0]*var[0])*(var[0]*var[0]))-1.6704000000000001e+00"; + } + + std::string + sym() const override + { + return "-1.6704+(1.3687)*c^4+(5.1622)*c^2-(2.7375)*c^3-(4.776)*c"; + } + + std::string + latex() const override + { + return "-1.6704+{(5.1622)} c^{2}-{(4.776)} c-{(2.7375)} c^{3}+{(1.3687)} c^{4}"; + } + + pfunct_faV_f * + clone() const override + { + return new pfunct_faV_f(*this); + } + }; + + template + class pfunct_faV_grad_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const override + { + return 5.4748000000000001e+00 * (var[0] * var[0] * var[0]) + + -8.2125000000000004e+00 * (var[0] * var[0]) + + 1.0324400000000001e+01 * var[0] - 4.7759999999999998e+00; + } + + public: + pfunct_faV_grad_0() + { + this->_name = "pfunct_faV_grad_0"; + } + + std::string + csrc() const override + { + return " 5.4748000000000001e+00*(var[0]*var[0]*var[0])+-8.2125000000000004e+00*(" + "var[0]*var[0])+1.0324400000000001e+01*var[0]-4.7759999999999998e+00"; + } + + std::string + sym() const override + { + return "-4.776+(10.3244)*c+(5.4748)*c^3-(8.2125)*c^2"; + } + + std::string + latex() const override + { + return "-4.776+{(5.4748)} c^{3}-{(8.2125)} c^{2}+{(10.3244)} c"; + } + + pfunct_faV_grad_0 * + clone() const override + { + return new pfunct_faV_grad_0(*this); + } + }; + + template + class pfunct_faV_hess_0_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const override + { + return 1.6424399999999999e+01 * (var[0] * var[0]) + + -1.6425000000000001e+01 * var[0] + 1.0324400000000001e+01; + } + + public: + pfunct_faV_hess_0_0() + { + this->_name = "pfunct_faV_hess_0_0"; + } + + std::string + csrc() const override + { + return " 1.6424399999999999e+01*(var[0]*var[0])+-1.6425000000000001e+01*var[0]+1." + "0324400000000001e+01"; + } + + std::string + sym() const override + { + return "10.3244-(16.425)*c+(16.4244)*c^2"; + } + + std::string + latex() const override + { + return "10.3244+{(16.4244)} c^{2}-{(16.425)} c"; + } + + pfunct_faV_hess_0_0 * + clone() const override + { + return new pfunct_faV_hess_0_0(*this); + } + }; + + template + class pfunct_faV : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_faV() + { + construct(); + } + + pfunct_faV(const pfunct_faV &RHS) + : PFuncBase(RHS) + { + construct(false); + + _val = RHS._val->clone(); + _grad_val[0] = RHS._grad_val[0]->clone(); + _hess_val[0][0] = RHS._hess_val[0][0]->clone(); + } + + pfunct_faV & + operator=(pfunct_faV RHS) + { + using std::swap; + + swap(_val, RHS._val); + swap(_grad_val[0], RHS._grad_val[0]); + swap(_hess_val[0][0], RHS._hess_val[0][0]); + + return *this; + } + + ~pfunct_faV() + { + delete _val; + + delete _grad_val[0]; + delete[] _grad_val; + + delete _hess_val[0][0]; + delete[] _hess_val[0]; + delete[] _hess_val; + } + + pfunct_faV * + clone() const override + { + return new pfunct_faV(*this); + } + + PSimpleFunction + simplefunction() const override + { + return PSimpleFunction(*_val); + } + + PSimpleFunction + grad_simplefunction(size_type di) const override + { + return PSimpleFunction(*_grad_val[di]); + } + + PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const override + { + return PSimpleFunction(*_hess_val[di][dj]); + } + + double + operator()(const VarContainer &var) override + { + return (*_val)(var); + } + + double + grad(const VarContainer &var, size_type di) override + { + return (*_grad_val[di])(var); + } + + double + hess(const VarContainer &var, size_type di, size_type dj) override + { + return (*_hess_val[di][dj])(var); + } + + void + eval(const VarContainer &var) override + { + (*_val)(var); + } + + void + eval_grad(const VarContainer &var) override + { + (*_grad_val[0])(var); + } + + void + eval_hess(const VarContainer &var) override + { + (*_hess_val[0][0])(var); + } + + double + operator()() const override + { + return (*_val)(); + } + + double + grad(size_type di) const override + { + return (*_grad_val[di])(); + } + + double + hess(size_type di, size_type dj) const override + { + return (*_hess_val[di][dj])(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_faV"; + this->_var_name.clear(); + this->_var_name.push_back("c"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + _grad_val = new PSimpleBase *[1]; + + _hess_val = new PSimpleBase **[1]; + _hess_val[0] = new PSimpleBase *[1]; + + if (!allocate) + return; + + _val = new pfunct_faV_f(); + + _grad_val[0] = new pfunct_faV_grad_0(); + + _hess_val[0][0] = new pfunct_faV_hess_0_0(); + } + }; + +} // namespace PRISMS +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh new file mode 100644 index 000000000..dcf1f9ef2 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/PLibrary/pfunct_fbV.hh @@ -0,0 +1,291 @@ +// created: 2017-1-11 12:28:13 +// version: master +// url: https://github.com/bpuchala/IntegrationToolsWriter.git +// commit: 13e063c3ac8e8911a726a243fdbd68f291cc58cc + +#ifndef pfunct_fbV_HH +#define pfunct_fbV_HH + +#include +#include +#include + +namespace PRISMS +{ + template + class pfunct_fbV_f : public PSimpleBase + { + double + eval(const VarContainer &var) const override + { + return -5.9745999999999997e+00 * var[0] + + 5.0000000000000000e+00 * (var[0] * var[0]) - 1.5924000000000000e+00; + } + + public: + pfunct_fbV_f() + { + this->_name = "pfunct_fbV_f"; + } + + std::string + csrc() const override + { + return " -5.9745999999999997e+00*var[0]+5.0000000000000000e+00*(var[0]*var[0])-1." + "5924000000000000e+00"; + } + + std::string + sym() const override + { + return "-1.5924-(5.9746)*c+(5.0)*c^2"; + } + + std::string + latex() const override + { + return "-1.5924+{(5.0)} c^{2}-{(5.9746)} c"; + } + + pfunct_fbV_f * + clone() const override + { + return new pfunct_fbV_f(*this); + } + }; + + template + class pfunct_fbV_grad_0 : public PSimpleBase + { + double + eval(const VarContainer &var) const override + { + return 1.0000000000000000e+01 * var[0] - 5.9745999999999997e+00; + } + + public: + pfunct_fbV_grad_0() + { + this->_name = "pfunct_fbV_grad_0"; + } + + std::string + csrc() const override + { + return " 1.0000000000000000e+01*var[0]-5.9745999999999997e+00"; + } + + std::string + sym() const override + { + return "-5.9746+(10.0)*c"; + } + + std::string + latex() const override + { + return "-5.9746+{(10.0)} c"; + } + + pfunct_fbV_grad_0 * + clone() const override + { + return new pfunct_fbV_grad_0(*this); + } + }; + + template + class pfunct_fbV_hess_0_0 : public PSimpleBase + { + double + eval([[maybe_unused]] const VarContainer &var) const override + { + return 1.0000000000000000e+01; + } + + public: + pfunct_fbV_hess_0_0() + { + this->_name = "pfunct_fbV_hess_0_0"; + } + + std::string + csrc() const override + { + return "1.0000000000000000e+01"; + } + + std::string + sym() const override + { + return "10.0"; + } + + std::string + latex() const override + { + return "10.0"; + } + + pfunct_fbV_hess_0_0 * + clone() const override + { + return new pfunct_fbV_hess_0_0(*this); + } + }; + + template + class pfunct_fbV : public PFuncBase + { + public: + typedef typename PFuncBase::size_type size_type; + + PSimpleBase *_val; + PSimpleBase **_grad_val; + PSimpleBase ***_hess_val; + + pfunct_fbV() + { + construct(); + } + + pfunct_fbV(const pfunct_fbV &RHS) + : PFuncBase(RHS) + { + construct(false); + + _val = RHS._val->clone(); + _grad_val[0] = RHS._grad_val[0]->clone(); + _hess_val[0][0] = RHS._hess_val[0][0]->clone(); + } + + pfunct_fbV & + operator=(pfunct_fbV RHS) + { + using std::swap; + + swap(_val, RHS._val); + swap(_grad_val[0], RHS._grad_val[0]); + swap(_hess_val[0][0], RHS._hess_val[0][0]); + + return *this; + } + + ~pfunct_fbV() + { + delete _val; + + delete _grad_val[0]; + delete[] _grad_val; + + delete _hess_val[0][0]; + delete[] _hess_val[0]; + delete[] _hess_val; + } + + pfunct_fbV * + clone() const override + { + return new pfunct_fbV(*this); + } + + PSimpleFunction + simplefunction() const override + { + return PSimpleFunction(*_val); + } + + PSimpleFunction + grad_simplefunction(size_type di) const override + { + return PSimpleFunction(*_grad_val[di]); + } + + PSimpleFunction + hess_simplefunction(size_type di, size_type dj) const override + { + return PSimpleFunction(*_hess_val[di][dj]); + } + + double + operator()(const VarContainer &var) override + { + return (*_val)(var); + } + + double + grad(const VarContainer &var, size_type di) override + { + return (*_grad_val[di])(var); + } + + double + hess(const VarContainer &var, size_type di, size_type dj) override + { + return (*_hess_val[di][dj])(var); + } + + void + eval(const VarContainer &var) override + { + (*_val)(var); + } + + void + eval_grad(const VarContainer &var) override + { + (*_grad_val[0])(var); + } + + void + eval_hess(const VarContainer &var) override + { + (*_hess_val[0][0])(var); + } + + double + operator()() const override + { + return (*_val)(); + } + + double + grad(size_type di) const override + { + return (*_grad_val[di])(); + } + + double + hess(size_type di, size_type dj) const override + { + return (*_hess_val[di][dj])(); + } + + private: + void + construct(bool allocate = true) + { + this->_name = "pfunct_fbV"; + this->_var_name.clear(); + this->_var_name.push_back("c"); + this->_var_description.clear(); + this->_var_description.push_back("concentration"); + + _grad_val = new PSimpleBase *[1]; + + _hess_val = new PSimpleBase **[1]; + _hess_val[0] = new PSimpleBase *[1]; + + if (!allocate) + return; + + _val = new pfunct_fbV_f(); + + _grad_val[0] = new pfunct_fbV_grad_0(); + + _hess_val[0][0] = new pfunct_fbV_hess_0_0(); + } + }; + +} // namespace PRISMS +#endif diff --git a/automatic_tests/precipitateEvolution_pfunction/customPDE.h b/automatic_tests/precipitateEvolution_pfunction/customPDE.h new file mode 100644 index 000000000..2985af610 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/customPDE.h @@ -0,0 +1,129 @@ +#include + +using namespace dealii; + +// Header files for PFunctions +typedef VectorizedArray scalarvalueType; +#include "PLibrary/PLibrary.cc" +#include "PLibrary/PLibrary.hh" + +#include + +// Declare the PFunctions to be used +PFunctions::pFunction pfunct_McV("pfunct_McV"), pfunct_Mn1V("pfunct_Mn1V"), + pfunct_Mn2V("pfunct_Mn1V"), pfunct_Mn3V("pfunct_Mn1V"), pfunct_faV("pfunct_faV"), + pfunct_fbV("pfunct_fbV"); + +template +class customPDE : public MatrixFreePDE +{ +public: + customPDE(userInputParameters _userInputs) + : MatrixFreePDE(_userInputs) + , userInputs(_userInputs) {}; + // Function to set the initial conditions (in ICs_and_BCs.h) + void + setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) override; + + // Function to set the non-uniform Dirichlet boundary conditions (in + // ICs_and_BCs.h) + void + setNonUniformDirichletBCs([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) override; + +private: +#include + + const userInputParameters userInputs; + + // Function to set the RHS of the governing equations for explicit time + // dependent equations (in equations.h) + void + explicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the RHS of the governing equations for all other equations + // (in equations.h) + void + nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the LHS of the governing equations (in equations.h) + void + equationLHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + +// Function to set postprocessing expressions (in postprocess.h) +#ifdef POSTPROCESS_FILE_EXISTS + void + postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; +#endif + +// Function to set the nucleation probability (in nucleation.h) +#ifdef NUCLEATION_FILE_EXISTS + double + getNucleationProbability([[maybe_unused]] variableValueContainer variable_value, + [[maybe_unused]] double dV) const override; +#endif + + // ================================================================ + // Methods specific to this subclass + // ================================================================ + + // ================================================================ + // Model constants specific to this subclass + // ================================================================ + + Tensor<2, dim> Kn1 = userInputs.get_model_constant_rank_2_tensor("Kn1"); + Tensor<2, dim> Kn2 = userInputs.get_model_constant_rank_2_tensor("Kn2"); + Tensor<2, dim> Kn3 = userInputs.get_model_constant_rank_2_tensor("Kn3"); + bool n_dependent_stiffness = + userInputs.get_model_constant_bool("n_dependent_stiffness"); + Tensor<2, dim> sfts_linear1 = + userInputs.get_model_constant_rank_2_tensor("sfts_linear1"); + Tensor<2, dim> sfts_const1 = userInputs.get_model_constant_rank_2_tensor("sfts_const1"); + Tensor<2, dim> sfts_linear2 = + userInputs.get_model_constant_rank_2_tensor("sfts_linear2"); + Tensor<2, dim> sfts_const2 = userInputs.get_model_constant_rank_2_tensor("sfts_const2"); + Tensor<2, dim> sfts_linear3 = + userInputs.get_model_constant_rank_2_tensor("sfts_linear3"); + Tensor<2, dim> sfts_const3 = userInputs.get_model_constant_rank_2_tensor("sfts_const3"); + double A4 = userInputs.get_model_constant_double("A4"); + double A3 = userInputs.get_model_constant_double("A3"); + double A2 = userInputs.get_model_constant_double("A2"); + double A1 = userInputs.get_model_constant_double("A1"); + double A0 = userInputs.get_model_constant_double("A0"); + double B2 = userInputs.get_model_constant_double("B2"); + double B1 = userInputs.get_model_constant_double("B1"); + double B0 = userInputs.get_model_constant_double("B0"); + + const static unsigned int CIJ_tensor_size = 2 * dim - 1 + dim / 3; + Tensor<2, CIJ_tensor_size> CIJ_Mg = + userInputs.get_model_constant_elasticity_tensor("CIJ_Mg"); + Tensor<2, CIJ_tensor_size> CIJ_Beta = + userInputs.get_model_constant_elasticity_tensor("CIJ_Beta"); + + // ================================================================ +}; diff --git a/automatic_tests/precipitateEvolution_pfunction/equations.cc b/automatic_tests/precipitateEvolution_pfunction/equations.cc new file mode 100644 index 000000000..edaa6603b --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/equations.cc @@ -0,0 +1,558 @@ +// ================================================================================= +// Set the attributes of the primary field variables +// ================================================================================= +// This function sets attributes for each variable/equation in the app. The +// attributes are set via standardized function calls. The first parameter for +// each function call is the variable index (starting at zero). The first set of +// variable/equation attributes are the variable name (any string), the variable +// type (SCALAR/VECTOR), and the equation type (EXPLICIT_TIME_DEPENDENT/ +// TIME_INDEPENDENT/AUXILIARY). The next set of attributes describe the +// dependencies for the governing equation on the values and derivatives of the +// other variables for the value term and gradient term of the RHS and the LHS. +// The final pair of attributes determine whether a variable represents a field +// that can nucleate and whether the value of the field is needed for nucleation +// rate calculations. + +void +variableAttributeLoader::loadVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "c"); + set_variable_type(0, SCALAR); + set_variable_equation_type(0, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(0, "c"); + set_dependencies_gradient_term_RHS( + 0, + "c, grad(c), n1, grad(n1), n2, grad(n2), n3, grad(n3), grad(u), hess(u)"); + + // Variable 1 + set_variable_name(1, "n1"); + set_variable_type(1, SCALAR); + set_variable_equation_type(1, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(1, "c, n1, n2, n3, grad(u)"); + set_dependencies_gradient_term_RHS(1, "grad(n1)"); + + // Variable 2 + set_variable_name(2, "n2"); + set_variable_type(2, SCALAR); + set_variable_equation_type(2, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(2, "c, n1, n2, n3, grad(u)"); + set_dependencies_gradient_term_RHS(2, "grad(n2)"); + + // Variable 3 + set_variable_name(3, "n3"); + set_variable_type(3, SCALAR); + set_variable_equation_type(3, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(3, "c, n1, n2, n3, grad(u)"); + set_dependencies_gradient_term_RHS(3, "grad(n3)"); + + // Variable 2 + set_variable_name(4, "u"); + set_variable_type(4, VECTOR); + set_variable_equation_type(4, TIME_INDEPENDENT); + + set_dependencies_value_term_RHS(4, ""); + set_dependencies_gradient_term_RHS(4, "c, n1, n2, n3, grad(u)"); + set_dependencies_value_term_LHS(4, ""); + set_dependencies_gradient_term_LHS(4, "n1, n2, n3, grad(change(u))"); +} + +// ============================================================================================= +// explicitEquationRHS (needed only if one or more equation is explict time +// dependent) +// ============================================================================================= +// This function calculates the right-hand-side of the explicit time-dependent +// equations for each variable. It takes "variable_list" as an input, which is a +// list of the value and derivatives of each of the variables at a specific +// quadrature point. The (x,y,z) location of that quadrature point is given by +// "q_point_loc". The function outputs two terms to variable_list -- one +// proportional to the test function and one proportional to the gradient of the +// test function. The index for each variable in this list corresponds to the +// index given at the top of this file. + +template +void +customPDE::explicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The concentration and its derivatives + scalarvalueType c = variable_list.get_scalar_value(0); + scalargradType cx = variable_list.get_scalar_gradient(0); + + // The first order parameter and its derivatives + scalarvalueType n1 = variable_list.get_scalar_value(1); + scalargradType n1x = variable_list.get_scalar_gradient(1); + + // The second order parameter and its derivatives + scalarvalueType n2 = variable_list.get_scalar_value(2); + scalargradType n2x = variable_list.get_scalar_gradient(2); + + // The third order parameter and its derivatives + scalarvalueType n3 = variable_list.get_scalar_value(3); + scalargradType n3x = variable_list.get_scalar_gradient(3); + + // The derivative of the displacement vector + vectorgradType ux = variable_list.get_vector_gradient(4); + + // --- Setting the expressions for the terms in the governing equations --- + + vectorhessType uxx; + + bool c_dependent_misfit = false; + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + if (std::abs(sfts_linear1[i][j]) > 1.0e-12) + { + c_dependent_misfit = true; + } + } + } + + if (c_dependent_misfit == true) + { + uxx = variable_list.get_vector_hessian(4); + } + + // Cahn-Hilliard mobility + scalarvalueType McV = pfunct_McV.val(c); + + // Allen-Cahn mobilities + scalarvalueType Mn1V = pfunct_Mn1V.val(c); + scalarvalueType Mn2V = pfunct_Mn2V.val(c); + scalarvalueType Mn3V = pfunct_Mn3V.val(c); + + // Free energy expressions and interpolation functions + scalarvalueType faV = pfunct_faV.val(c); + scalarvalueType facV = pfunct_faV.grad(c, 0); + scalarvalueType faccV = pfunct_faV.hess(c, 0, 0); + scalarvalueType fbV = pfunct_fbV.val(c); + scalarvalueType fbcV = pfunct_fbV.grad(c, 0); + scalarvalueType fbccV = pfunct_fbV.hess(c, 0, 0); + + scalarvalueType h1V = + (10.0 * n1 * n1 * n1 - 15.0 * n1 * n1 * n1 * n1 + 6.0 * n1 * n1 * n1 * n1 * n1); + scalarvalueType h2V = + (10.0 * n2 * n2 * n2 - 15.0 * n2 * n2 * n2 * n2 + 6.0 * n2 * n2 * n2 * n2 * n2); + scalarvalueType h3V = + (10.0 * n3 * n3 * n3 - 15.0 * n3 * n3 * n3 * n3 + 6.0 * n3 * n3 * n3 * n3 * n3); + scalarvalueType hn1V = + (30.0 * n1 * n1 - 60.0 * n1 * n1 * n1 + 30.0 * n1 * n1 * n1 * n1); + scalarvalueType hn2V = + (30.0 * n2 * n2 - 60.0 * n2 * n2 * n2 + 30.0 * n2 * n2 * n2 * n2); + scalarvalueType hn3V = + (30.0 * n3 * n3 - 60.0 * n3 * n3 * n3 + 30.0 * n3 * n3 * n3 * n3); + + // Calculate the stress-free transformation strain and its derivatives at the + // quadrature point + Tensor<2, dim, VectorizedArray> sfts1, sfts1c, sfts1cc, sfts2, sfts2c, sfts2cc, + sfts3, sfts3c, sfts3cc; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts1[i][j] = constV(sfts_linear1[i][j]) * c + constV(sfts_const1[i][j]); + sfts1c[i][j] = constV(sfts_linear1[i][j]); + sfts1cc[i][j] = constV(0.0); + + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts2[i][j] = constV(sfts_linear2[i][j]) * c + constV(sfts_const2[i][j]); + sfts2c[i][j] = constV(sfts_linear1[i][j]); + sfts2cc[i][j] = constV(0.0); + + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts3[i][j] = constV(sfts_linear3[i][j]) * c + constV(sfts_const3[i][j]); + sfts3c[i][j] = constV(sfts_linear3[i][j]); + sfts3cc[i][j] = constV(0.0); + } + } + + // compute E2=(E-E0) + VectorizedArray E2[dim][dim], S[dim][dim]; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + E2[i][j] = constV(0.5) * (ux[i][j] + ux[j][i]) - + (sfts1[i][j] * h1V + sfts2[i][j] * h2V + sfts3[i][j] * h3V); + } + } + + // compute stress + // S=C*(E-E0) + // Compute stress tensor (which is equal to the residual, Rux) + VectorizedArray CIJ_combined[CIJ_tensor_size][CIJ_tensor_size]; + + if (n_dependent_stiffness == true) + { + VectorizedArray sum_hV; + sum_hV = h1V + h2V + h3V; + for (unsigned int i = 0; i < 2 * dim - 1 + dim / 3; i++) + { + for (unsigned int j = 0; j < 2 * dim - 1 + dim / 3; j++) + { + CIJ_combined[i][j] = + CIJ_Mg[i][j] * (constV(1.0) - sum_hV) + CIJ_Beta[i][j] * sum_hV; + } + } + computeStress(CIJ_combined, E2, S); + } + else + { + computeStress(CIJ_Mg, E2, S); + } + + // Compute one of the stress terms in the order parameter chemical potential, + // nDependentMisfitACp = C*(E-E0)*(E0_p*Hn) + VectorizedArray nDependentMisfitAC1 = constV(0.0); + VectorizedArray nDependentMisfitAC2 = constV(0.0); + VectorizedArray nDependentMisfitAC3 = constV(0.0); + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + nDependentMisfitAC1 += S[i][j] * (sfts1[i][j]); + nDependentMisfitAC2 += S[i][j] * (sfts2[i][j]); + nDependentMisfitAC3 += S[i][j] * (sfts3[i][j]); + } + } + + nDependentMisfitAC1 *= -hn1V; + nDependentMisfitAC2 *= -hn2V; + nDependentMisfitAC3 *= -hn3V; + + // Compute the other stress term in the order parameter chemical potential, + // heterMechACp = 0.5*Hn*(C_beta-C_alpha)*(E-E0)*(E-E0) + VectorizedArray heterMechAC1 = constV(0.0); + VectorizedArray heterMechAC2 = constV(0.0); + VectorizedArray heterMechAC3 = constV(0.0); + VectorizedArray S2[dim][dim]; + + if (n_dependent_stiffness == true) + { + // computeStress(CIJ_diff, E2, S2); + computeStress(CIJ_Beta - CIJ_Mg, E2, S2); + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + heterMechAC1 += S2[i][j] * E2[i][j]; + } + } + // Aside from HnpV, heterMechAC1, heterMechAC2, and heterMechAC3 are equal + heterMechAC2 = 0.5 * hn2V * heterMechAC1; + heterMechAC3 = 0.5 * hn3V * heterMechAC1; + + heterMechAC1 = 0.5 * hn1V * heterMechAC1; + } + + // compute the stress term in the gradient of the concentration chemical + // potential, grad_mu_el = [C*(E-E0)*E0c]x, must be a vector with length dim + scalargradType grad_mu_el; + + if (c_dependent_misfit == true) + { + VectorizedArray E3[dim][dim], S3[dim][dim]; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + E3[i][j] = -(sfts1c[i][j] * h1V + sfts2c[i][j] * h2V + sfts3c[i][j] * h3V); + } + } + + if (n_dependent_stiffness == true) + { + computeStress(CIJ_combined, E3, S3); + } + else + { + computeStress(CIJ_Mg, E3, S3); + } + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + for (unsigned int k = 0; k < dim; k++) + { + grad_mu_el[k] += + S3[i][j] * + (constV(0.5) * (uxx[i][j][k] + uxx[j][i][k]) + E3[i][j] * cx[k] - + (sfts1[i][j] * hn1V * n1x[k] + sfts2[i][j] * hn2V * n2x[k] + + sfts3[i][j] * hn3V * n3x[k])); + + grad_mu_el[k] += + -S[i][j] * + (sfts1c[i][j] * hn1V * n1x[k] + sfts2c[i][j] * hn2V * n2x[k] + + sfts3c[i][j] * hn3V * n3x[k] + + (sfts1cc[i][j] * h1V + sfts2cc[i][j] * h2V + sfts3cc[i][j] * h3V) * + cx[k]); + + if (n_dependent_stiffness == true) + { + grad_mu_el[k] += S2[i][j] * E3[i][j] * + (hn1V * n1x[k] + hn2V * n2x[k] + hn3V * n3x[k]); + } + } + } + } + } + + // compute K*nx + scalargradType Knx1, Knx2, Knx3; + for (unsigned int a = 0; a < dim; a++) + { + Knx1[a] = 0.0; + Knx2[a] = 0.0; + Knx3[a] = 0.0; + for (unsigned int b = 0; b < dim; b++) + { + Knx1[a] += constV(Kn1[a][b]) * n1x[b]; + Knx2[a] += constV(Kn2[a][b]) * n2x[b]; + Knx3[a] += constV(Kn3[a][b]) * n3x[b]; + } + } + + // The terms in the govering equations + scalarvalueType eq_c = (c); + scalargradType eqx_c_temp = + (cx * ((1.0 - h1V - h2V - h3V) * faccV + (h1V + h2V + h3V) * fbccV) + + n1x * ((fbcV - facV) * hn1V) + n2x * ((fbcV - facV) * hn2V) + + n3x * ((fbcV - facV) * hn3V) + grad_mu_el); + scalargradType eqx_c = (constV(-userInputs.dtValue) * McV * eqx_c_temp); + + scalarvalueType eq_n1 = + (n1 - constV(userInputs.dtValue) * Mn1V * + ((fbV - faV) * hn1V + nDependentMisfitAC1 + heterMechAC1)); + scalarvalueType eq_n2 = + (n2 - constV(userInputs.dtValue) * Mn2V * + ((fbV - faV) * hn2V + nDependentMisfitAC2 + heterMechAC2)); + scalarvalueType eq_n3 = + (n3 - constV(userInputs.dtValue) * Mn3V * + ((fbV - faV) * hn3V + nDependentMisfitAC3 + heterMechAC3)); + scalargradType eqx_n1 = (constV(-userInputs.dtValue) * Mn1V * Knx1); + scalargradType eqx_n2 = (constV(-userInputs.dtValue) * Mn2V * Knx2); + scalargradType eqx_n3 = (constV(-userInputs.dtValue) * Mn3V * Knx3); + + // --- Submitting the terms for the governing equations --- + + variable_list.set_scalar_value_term_RHS(0, eq_c); + variable_list.set_scalar_gradient_term_RHS(0, eqx_c); + + variable_list.set_scalar_value_term_RHS(1, eq_n1); + variable_list.set_scalar_gradient_term_RHS(1, eqx_n1); + + variable_list.set_scalar_value_term_RHS(2, eq_n2); + variable_list.set_scalar_gradient_term_RHS(2, eqx_n2); + + variable_list.set_scalar_value_term_RHS(3, eq_n3); + variable_list.set_scalar_gradient_term_RHS(3, eqx_n3); +} + +// ============================================================================================= +// nonExplicitEquationRHS (needed only if one or more equation is time +// independent or auxiliary) +// ============================================================================================= +// This function calculates the right-hand-side of all of the equations that are +// not explicit time-dependent equations. It takes "variable_list" as an input, +// which is a list of the value and derivatives of each of the variables at a +// specific quadrature point. The (x,y,z) location of that quadrature point is +// given by "q_point_loc". The function outputs two terms to variable_list -- +// one proportional to the test function and one proportional to the gradient of +// the test function. The index for each variable in this list corresponds to +// the index given at the top of this file. + +template +void +customPDE::nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The concentration and its derivatives + scalarvalueType c = variable_list.get_scalar_value(0); + + // The first order parameter and its derivatives + scalarvalueType n1 = variable_list.get_scalar_value(1); + + // The second order parameter and its derivatives + scalarvalueType n2 = variable_list.get_scalar_value(2); + + // The third order parameter and its derivatives + scalarvalueType n3 = variable_list.get_scalar_value(3); + + // The derivative of the displacement vector + vectorgradType ux = variable_list.get_vector_gradient(4); + + // --- Setting the expressions for the terms in the governing equations --- + + // Interpolation functions + scalarvalueType h1V = + (10.0 * n1 * n1 * n1 - 15.0 * n1 * n1 * n1 * n1 + 6.0 * n1 * n1 * n1 * n1 * n1); + scalarvalueType h2V = + (10.0 * n2 * n2 * n2 - 15.0 * n2 * n2 * n2 * n2 + 6.0 * n2 * n2 * n2 * n2 * n2); + scalarvalueType h3V = + (10.0 * n3 * n3 * n3 - 15.0 * n3 * n3 * n3 * n3 + 6.0 * n3 * n3 * n3 * n3 * n3); + + // Calculate the stress-free transformation strain and its derivatives at the + // quadrature point + Tensor<2, dim, VectorizedArray> sfts1, sfts2, sfts3; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts1[i][j] = constV(sfts_linear1[i][j]) * c + constV(sfts_const1[i][j]); + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts2[i][j] = constV(sfts_linear2[i][j]) * c + constV(sfts_const2[i][j]); + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts3[i][j] = constV(sfts_linear3[i][j]) * c + constV(sfts_const3[i][j]); + } + } + + // compute E2=(E-E0) + VectorizedArray E2[dim][dim], S[dim][dim]; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + E2[i][j] = constV(0.5) * (ux[i][j] + ux[j][i]) - + (sfts1[i][j] * h1V + sfts2[i][j] * h2V + sfts3[i][j] * h3V); + } + } + + // compute stress + // S=C*(E-E0) + // Compute stress tensor (which is equal to the residual, Rux) + VectorizedArray CIJ_combined[CIJ_tensor_size][CIJ_tensor_size]; + + if (n_dependent_stiffness == true) + { + VectorizedArray sum_hV; + sum_hV = h1V + h2V + h3V; + for (unsigned int i = 0; i < 2 * dim - 1 + dim / 3; i++) + { + for (unsigned int j = 0; j < 2 * dim - 1 + dim / 3; j++) + { + CIJ_combined[i][j] = + CIJ_Mg[i][j] * (constV(1.0) - sum_hV) + CIJ_Beta[i][j] * sum_hV; + } + } + computeStress(CIJ_combined, E2, S); + } + else + { + computeStress(CIJ_Mg, E2, S); + } + + vectorgradType eqx_u; + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + eqx_u[i][j] = -S[i][j]; + } + } + + // --- Submitting the terms for the governing equations --- + + variable_list.set_vector_gradient_term_RHS(4, eqx_u); +} + +// ============================================================================================= +// equationLHS (needed only if at least one equation is time independent) +// ============================================================================================= +// This function calculates the left-hand-side of time-independent equations. It +// takes "variable_list" as an input, which is a list of the value and +// derivatives of each of the variables at a specific quadrature point. The +// (x,y,z) location of that quadrature point is given by "q_point_loc". The +// function outputs two terms to variable_list -- one proportional to the test +// function and one proportional to the gradient of the test function -- for the +// left-hand-side of the equation. The index for each variable in this list +// corresponds to the index given at the top of this file. If there are multiple +// elliptic equations, conditional statements should be sed to ensure that the +// correct residual is being submitted. The index of the field being solved can +// be accessed by "this->currentFieldIndex". + +template +void +customPDE::equationLHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // n1 + scalarvalueType n1 = variable_list.get_scalar_value(1); + + // n2 + scalarvalueType n2 = variable_list.get_scalar_value(2); + + // n3 + scalarvalueType n3 = variable_list.get_scalar_value(3); + + // u + vectorgradType Dux = variable_list.get_change_in_vector_gradient(4); + + // --- Setting the expressions for the terms in the governing equations --- + + vectorgradType eqx_Du; + + // Interpolation functions + + scalarvalueType h1V = + (10.0 * n1 * n1 * n1 - 15.0 * n1 * n1 * n1 * n1 + 6.0 * n1 * n1 * n1 * n1 * n1); + scalarvalueType h2V = + (10.0 * n2 * n2 * n2 - 15.0 * n2 * n2 * n2 * n2 + 6.0 * n2 * n2 * n2 * n2 * n2); + scalarvalueType h3V = + (10.0 * n3 * n3 * n3 - 15.0 * n3 * n3 * n3 * n3 + 6.0 * n3 * n3 * n3 * n3 * n3); + + // Take advantage of E being simply 0.5*(ux + transpose(ux)) and use the + // dealii "symmetrize" function + Tensor<2, dim, VectorizedArray> E; + E = symmetrize(Dux); + + // Compute stress tensor (which is equal to the residual, Rux) + if (n_dependent_stiffness == true) + { + Tensor<2, CIJ_tensor_size, VectorizedArray> CIJ_combined; + CIJ_combined = + CIJ_Mg * (constV(1.0) - h1V - h2V - h3V) + CIJ_Beta * (h1V + h2V + h3V); + + computeStress(CIJ_combined, E, eqx_Du); + } + else + { + computeStress(CIJ_Mg, E, eqx_Du); + } + + // --- Submitting the terms for the governing equations --- + + variable_list.set_vector_gradient_term_LHS(4, eqx_Du); +} diff --git a/automatic_tests/precipitateEvolution_pfunction/gold_integratedFields.txt b/automatic_tests/precipitateEvolution_pfunction/gold_integratedFields.txt new file mode 100644 index 000000000..a540bf6be --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/gold_integratedFields.txt @@ -0,0 +1,7 @@ +0 f_tot -2961.102683 +0.6 f_tot -2965.562629 +1.2 f_tot -2965.677538 +1.8 f_tot -2965.781569 +2.4 f_tot -2965.884576 +3 f_tot -2965.976361 +3.6 f_tot -2966.061707 diff --git a/automatic_tests/precipitateEvolution_pfunction/parameters.prm b/automatic_tests/precipitateEvolution_pfunction/parameters.prm new file mode 100644 index 000000000..b9b19c596 --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/parameters.prm @@ -0,0 +1,127 @@ +# ================================================================================= +# Set the number of dimensions (1, 2, or 3 for a 1D, 2D, or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# ================================================================================= +set Domain size X = 40.0 +set Domain size Y = 40.0 +set Domain size Z = 40.0 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 3 +set Subdivisions Y = 3 +set Subdivisions Z = 3 + +set Refine factor = 5 + +set Element degree = 2 + +# ================================================================================= +# Set the adaptive mesh refinement parameters +# ================================================================================= +set Mesh adaptivity = true + +set Max refinement level = 5 +set Min refinement level = 1 + +set Steps between remeshing operations = 1000 + +subsection Refinement criterion: n1 + set Criterion type = VALUE + set Value lower bound = 0.01 + set Value upper bound = 0.99 +end + +subsection Refinement criterion: n2 + set Criterion type = VALUE + set Value lower bound = 0.01 + set Value upper bound = 0.99 +end + +subsection Refinement criterion: n3 + set Criterion type = VALUE + set Value lower bound = 0.01 + set Value upper bound = 0.99 +end + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 6.0e-4 + +set Number of time steps = 6000 + +# ================================================================================= +# Set the linear solver parameters +# ================================================================================= + +subsection Linear solver parameters: u + set Tolerance type = ABSOLUTE_RESIDUAL + + set Tolerance value = 5e-3 + + set Maximum linear solver iterations = 1000 +end + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 6 + +set Skip print steps = 1000 + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable c = NATURAL +set Boundary condition for variable n1 = NATURAL +set Boundary condition for variable n2 = NATURAL +set Boundary condition for variable n3 = NATURAL +set Boundary condition for variable u, x component = DIRICHLET: 0.0 +set Boundary condition for variable u, y component = DIRICHLET: 0.0 +set Boundary condition for variable u, z component = DIRICHLET: 0.0 + +#set Boundary condition for variable u, x component = DIRICHLET: 0.0, DIRICHLET: 1.0, NATURAL, NATURAL +#set Boundary condition for variable u, y component = NATURAL +#set Boundary condition for variable u, z component = NATURAL + +# ================================================================================= +# Set the model constants +# ================================================================================= +# The gradient energy coefficients +set Model constant Kn1 = ((0.03,0,0),(0,0.007,0),(0,0,1.0)), tensor +set Model constant Kn2 = ((0.01275,-0.009959,0),(-0.009959,0.02425,0),(0,0,1.0)), tensor +set Model constant Kn3 = ((0.01275,0.009959,0),(0.009959,0.02425,0),(0,0,1.0)), tensor + +# n_dependent_stiffness +set Model constant n_dependent_stiffness = true, bool + +# The linear and constant coefficients of the stress-free transformation strains +set Model constant sfts_linear1 = ((0,0,0),(0,0,0),(0,0,0)), tensor +set Model constant sfts_const1 = ((0.0345,0,0),(0,0.0185,0),(0,0,-0.00270)), tensor +set Model constant sfts_linear2 = ((0,0,0),(0,0,0),(0,0,0)), tensor +set Model constant sfts_const2 = ((0.0225,-0.0069,0),(-0.0069,0.0305,0),(0,0,-0.00270)), tensor +set Model constant sfts_linear3 = ((0,0,0),(0,0,0),(0,0,0)), tensor +set Model constant sfts_const3 = ((0.0225, 0.0069,0),(0.0069,0.0305,0),(0,0,-0.00270)), tensor + +# A4, A3, A2, A1, and A0 Mg-Y matrix free energy parameters +set Model constant A4 = 1.3687, double +set Model constant A3 = -2.7375, double +set Model constant A2 = 5.1622, double +set Model constant A1 = -4.776, double +set Model constant A0 = -1.6704, double + +# B2, B1, and B0 Mg-Y matrix free energy parameters +set Model constant B2 = 5.0, double +set Model constant B1 = -5.9746, double +set Model constant B0 = -1.5924, double + +set Model constant CIJ_Mg = (40.0,0.3), isotropic elastic constants +set Model constant CIJ_Beta = (50.0,0.3), isotropic elastic constants diff --git a/automatic_tests/precipitateEvolution_pfunction/postprocess.cc b/automatic_tests/precipitateEvolution_pfunction/postprocess.cc new file mode 100644 index 000000000..077889a5c --- /dev/null +++ b/automatic_tests/precipitateEvolution_pfunction/postprocess.cc @@ -0,0 +1,187 @@ +// ============================================================================================= +// loadPostProcessorVariableAttributes: Set the attributes of the postprocessing +// variables +// ============================================================================================= +// This function is analogous to 'loadVariableAttributes' in 'equations.h', but +// for the postprocessing expressions. It sets the attributes for each +// postprocessing expression, including its name, whether it is a vector or +// scalar (only scalars are supported at present), its dependencies on other +// variables and their derivatives, and whether to calculate an integral of the +// postprocessed quantity over the entire domain. Note: this function is not a +// member of customPDE. + +void +variableAttributeLoader::loadPostProcessorVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "f_tot"); + set_variable_type(0, SCALAR); + + set_dependencies_value_term_RHS( + 0, + "c, grad(c), n1, grad(n1), n2, grad(n2), n3, grad(n3), grad(u)"); + set_dependencies_gradient_term_RHS(0, ""); + + set_output_integral(0, true); +} + +// ================================================================================= +// Define the expressions for the post-processed fields +// ================================================================================= + +template +void +customPDE::postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The concentration and its derivatives + scalarvalueType c = variable_list.get_scalar_value(0); + + // The first order parameter and its derivatives + scalarvalueType n1 = variable_list.get_scalar_value(1); + scalargradType n1x = variable_list.get_scalar_gradient(1); + + // The second order parameter and its derivatives + scalarvalueType n2 = variable_list.get_scalar_value(2); + scalargradType n2x = variable_list.get_scalar_gradient(2); + + // The third order parameter and its derivatives + scalarvalueType n3 = variable_list.get_scalar_value(3); + scalargradType n3x = variable_list.get_scalar_gradient(3); + + // The derivative of the displacement vector + vectorgradType ux = variable_list.get_vector_gradient(4); + + // --- Setting the expressions for the terms in the postprocessing expressions + // --- + + scalarvalueType f_tot = constV(0.0); + + // Free energy expressions and interpolation functions + scalarvalueType faV = pfunct_faV.val(c); + scalarvalueType fbV = pfunct_fbV.val(c); + scalarvalueType h1V = + (10.0 * n1 * n1 * n1 - 15.0 * n1 * n1 * n1 * n1 + 6.0 * n1 * n1 * n1 * n1 * n1); + scalarvalueType h2V = + (10.0 * n2 * n2 * n2 - 15.0 * n2 * n2 * n2 * n2 + 6.0 * n2 * n2 * n2 * n2 * n2); + scalarvalueType h3V = + (10.0 * n3 * n3 * n3 - 15.0 * n3 * n3 * n3 * n3 + 6.0 * n3 * n3 * n3 * n3 * n3); + + scalarvalueType f_chem = + (constV(1.0) - (h1V + h2V + h3V)) * faV + (h1V + h2V + h3V) * fbV; + + scalarvalueType f_grad = constV(0.0); + + for (int i = 0; i < dim; i++) + { + for (int j = 0; j < dim; j++) + { + f_grad += constV(0.5 * Kn1[i][j]) * n1x[i] * n1x[j]; + } + } + + for (int i = 0; i < dim; i++) + { + for (int j = 0; j < dim; j++) + { + f_grad += constV(0.5 * Kn2[i][j]) * n2x[i] * n2x[j]; + } + } + + for (int i = 0; i < dim; i++) + { + for (int j = 0; j < dim; j++) + { + f_grad += constV(0.5 * Kn3[i][j]) * n3x[i] * n3x[j]; + } + } + + // Calculate the stress-free transformation strain and its derivatives at the + // quadrature point + Tensor<2, dim, VectorizedArray> sfts1, sfts1c, sfts1cc, sfts2, sfts2c, sfts2cc, + sfts3, sfts3c, sfts3cc; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts1[i][j] = constV(sfts_linear1[i][j]) * c + constV(sfts_const1[i][j]); + sfts1c[i][j] = constV(sfts_linear1[i][j]); + sfts1cc[i][j] = constV(0.0); + + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts2[i][j] = constV(sfts_linear2[i][j]) * c + constV(sfts_const2[i][j]); + sfts2c[i][j] = constV(sfts_linear1[i][j]); + sfts2cc[i][j] = constV(0.0); + + // Polynomial fits for the stress-free transformation strains, of the + // form: sfts = a_p * c + b_p + sfts3[i][j] = constV(sfts_linear3[i][j]) * c + constV(sfts_const3[i][j]); + sfts3c[i][j] = constV(sfts_linear3[i][j]); + sfts3cc[i][j] = constV(0.0); + } + } + + // compute E2=(E-E0) + VectorizedArray E2[dim][dim], S[dim][dim]; + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + // E2[i][j]= constV(0.5)*(ux[i][j]+ux[j][i])-( sfts1[i][j]*h1V + + // sfts2[i][j]*h2V + sfts3[i][j]*h3V); + E2[i][j] = constV(0.5) * (ux[i][j] + ux[j][i]) - + (sfts1[i][j] * h1V + sfts2[i][j] * h2V + sfts3[i][j] * h3V); + } + } + + // compute stress + // S=C*(E-E0) + VectorizedArray CIJ_combined[CIJ_tensor_size][CIJ_tensor_size]; + + if (n_dependent_stiffness == true) + { + VectorizedArray sum_hV; + sum_hV = h1V + h2V + h3V; + for (unsigned int i = 0; i < 2 * dim - 1 + dim / 3; i++) + { + for (unsigned int j = 0; j < 2 * dim - 1 + dim / 3; j++) + { + CIJ_combined[i][j] = + CIJ_Mg[i][j] * (constV(1.0) - sum_hV) + CIJ_Beta[i][j] * sum_hV; + } + } + computeStress(CIJ_combined, E2, S); + } + else + { + computeStress(CIJ_Mg, E2, S); + } + + scalarvalueType f_el = constV(0.0); + + for (unsigned int i = 0; i < dim; i++) + { + for (unsigned int j = 0; j < dim; j++) + { + f_el += constV(0.5) * S[i][j] * E2[i][j]; + } + } + + f_tot = f_chem + f_grad + f_el; + + // --- Submitting the terms for the postprocessing expressions --- + + pp_variable_list.set_scalar_value_term_RHS(0, f_tot); +} diff --git a/tests/automatic_tests/run_automatic_tests.py b/automatic_tests/run_automatic_tests.py similarity index 96% rename from tests/automatic_tests/run_automatic_tests.py rename to automatic_tests/run_automatic_tests.py index 691b9500c..c1973a1af 100644 --- a/tests/automatic_tests/run_automatic_tests.py +++ b/automatic_tests/run_automatic_tests.py @@ -149,7 +149,7 @@ def run_regression_test(application, new_gold_standard, test_dir): rel_diff = abs( (float(gold_last_energy) - float(last_energy)) / float(gold_last_energy) ) - test_passed = rel_diff < 1.0e-5 + test_passed = rel_diff < 1.0e-4 # Determine test result test_result = ( @@ -220,10 +220,26 @@ def run_regression_tests_in_parallel( "allenCahn", "cahnHilliard", "CHAC_anisotropyRegularized", + "CHiMaD_benchmark6a", + "corrosion_microgalvanic", "coupledCahnHilliardAllenCahn", + "grainGrowth", "precipitateEvolution", + "precipitateEvolution_pfunction", + "spinodalDecomposition", +] +getNewGoldStandardList = [ + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, ] -getNewGoldStandardList = [False, False, False, False, False] # Grab cpu information architecture, cpu_model, cpu_cores, cpu_max_freq, cpu_min_freq, hypervisor = ( diff --git a/automatic_tests/spinodalDecomposition/CMakeLists.txt b/automatic_tests/spinodalDecomposition/CMakeLists.txt new file mode 100644 index 000000000..82b46ea0a --- /dev/null +++ b/automatic_tests/spinodalDecomposition/CMakeLists.txt @@ -0,0 +1,165 @@ +## +# CMake script for the PRISMS-PF applications +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(myapp) + +message(STATUS "") +message(STATUS "=========================================================") +message(STATUS "Configuring PRISMS-PF application") +message(STATUS "=========================================================") +message(STATUS "") + +# Setup the build type (debug, release, debugrelease) +if("${CMAKE_BUILD_TYPE}" STREQUAL "") + set(CMAKE_BUILD_TYPE "DebugRelease" + CACHE STRING + "Choose the type of build, options are: Debug, Release and DebugRelease." + FORCE) +endif() + +# Convert build type into the debug and release builds, which may or may +# not be built. +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) + message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") +else() + message(FATAL_ERROR + "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_DEBUG "ON") +else() + set(PRISMS_PF_BUILD_DEBUG "OFF") +endif() + +if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR + "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") + set(PRISMS_PF_BUILD_RELEASE "ON") +else() + set(PRISMS_PF_BUILD_RELEASE "OFF") +endif() + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Make and ninja build options +if(CMAKE_GENERATOR MATCHES "Ninja") + set(_make_command "$ ninja") +else() + set(_make_command " $ make") +endif() + +# Debug and release targets +if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") + add_custom_target(release + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to RELEASE mode..." + ) + + add_custom_target(debug + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG mode..." + ) + + add_custom_target(debugrelease + COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . + COMMAND ${CMAKE_COMMAND} -E echo "***" + COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" + COMMAND ${CMAKE_COMMAND} -E echo "***" + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + VERBATIM + COMMENT "switching to DEBUG/RELEASE mode..." + ) +endif() + +# Add distclean target to clean build +add_custom_target(distclean + COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean + COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles + COMMAND ${CMAKE_COMMAND} -E remove + CMakeCache.txt cmake_install.cmake Makefile + build.ninja rules.ninja .ninja_deps .ninja_log + COMMENT "distclean invoked" + ) + +# Add postprocess.cc and nucleation.cc if they exist +if(EXISTS "postprocess.cc") + add_definitions(-DPOSTPROCESS_FILE_EXISTS) +endif() +if(EXISTS "nucleation.cc") + add_definitions(-DNUCLEATION_FILE_EXISTS) +endif() + +# Set location of files +include_directories(${CMAKE_SOURCE_DIR}/../../include) +include_directories(${CMAKE_SOURCE_DIR}/../../src) +include_directories(${CMAKE_SOURCE_DIR}) + +# Set the location of the main.cc file +set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") + add_executable(main_debug ${TARGET_SRC}) + set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) + deal_ii_setup_target(main_debug DEBUG) + target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../libprisms-pf-debug.a) +endif() + +if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") + add_executable(main_release ${TARGET_SRC}) + set_property(TARGET main_release PROPERTY OUTPUT_NAME main) + deal_ii_setup_target(main_release RELEASE) + target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../libprisms-pf-release.a) +endif() \ No newline at end of file diff --git a/automatic_tests/spinodalDecomposition/ICs_and_BCs.cc b/automatic_tests/spinodalDecomposition/ICs_and_BCs.cc new file mode 100644 index 000000000..95046044b --- /dev/null +++ b/automatic_tests/spinodalDecomposition/ICs_and_BCs.cc @@ -0,0 +1,62 @@ +// =========================================================================== +// FUNCTION FOR INITIAL CONDITIONS +// =========================================================================== + +template +void +customPDE::setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) +{ + // --------------------------------------------------------------------- + // ENTER THE INITIAL CONDITIONS HERE + // --------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the initial condition for each variable + // according to its variable index + + // The initial condition is a set of overlapping circles/spheres defined + // by a hyperbolic tangent function. The center of each circle/sphere is + // given by "center" and its radius is given by "radius". + + if (index == 0) + { + scalar_IC = c0 + 2.0 * icamplitude * (dist(rng) - 0.5); + } + else + { + scalar_IC = 0.0; + } + + // --------------------------------------------------------------------- +} + +// =========================================================================== +// FUNCTION FOR NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS +// =========================================================================== + +template +void +customPDE::setNonUniformDirichletBCs( + [[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) +{ + // -------------------------------------------------------------------------- + // ENTER THE NON-UNIFORM DIRICHLET BOUNDARY CONDITIONS HERE + // -------------------------------------------------------------------------- + // Enter the function describing conditions for the fields at point "p". + // Use "if" statements to set the boundary condition for each variable + // according to its variable index. This function can be left blank if there + // are no non-uniform Dirichlet boundary conditions. For BCs that change in + // time, you can access the current time through the variable "time". The + // boundary index can be accessed via the variable "direction", which starts + // at zero and uses the same order as the BC specification in parameters.in + // (i.e. left = 0, right = 1, bottom = 2, top = 3, front = 4, back = 5). + + // ------------------------------------------------------------------------- +} diff --git a/automatic_tests/spinodalDecomposition/customPDE.h b/automatic_tests/spinodalDecomposition/customPDE.h new file mode 100644 index 000000000..34ed959e4 --- /dev/null +++ b/automatic_tests/spinodalDecomposition/customPDE.h @@ -0,0 +1,113 @@ +#include +#include + +using namespace dealii; + +template +class customPDE : public MatrixFreePDE +{ +public: + // Constructor + customPDE(userInputParameters _userInputs) + : MatrixFreePDE(_userInputs) + , userInputs(_userInputs) + { + // Defining seed + unsigned seed = 123456789; + // Initializing distribution + dist = distribution(0.0, 1.0); + // Initializing random variable + rng = engine(seed); + }; + + // Function to set the initial conditions (in ICs_and_BCs.h) + void + setInitialCondition([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] double &scalar_IC, + [[maybe_unused]] Vector &vector_IC) override; + + // Function to set the non-uniform Dirichlet boundary conditions (in + // ICs_and_BCs.h) + void + setNonUniformDirichletBCs([[maybe_unused]] const Point &p, + [[maybe_unused]] const unsigned int index, + [[maybe_unused]] const unsigned int direction, + [[maybe_unused]] const double time, + [[maybe_unused]] double &scalar_BC, + [[maybe_unused]] Vector &vector_BC) override; + + typedef std::mt19937_64 engine; + typedef std::uniform_real_distribution distribution; + +private: +#include + + const userInputParameters userInputs; + + // Function to set the RHS of the governing equations for explicit time + // dependent equations (in equations.h) + void + explicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the RHS of the governing equations for all other equations + // (in equations.h) + void + nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + + // Function to set the LHS of the governing equations (in equations.h) + void + equationLHS( + [[maybe_unused]] variableContainer> + &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; + +// Function to set postprocessing expressions (in postprocess.h) +#ifdef POSTPROCESS_FILE_EXISTS + void + postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const override; +#endif + +// Function to set the nucleation probability (in nucleation.h) +#ifdef NUCLEATION_FILE_EXISTS + double + getNucleationProbability([[maybe_unused]] variableValueContainer variable_value, + [[maybe_unused]] double dV) const override; +#endif + + // ================================================================ + // Methods specific to this subclass + // ================================================================ + + // ================================================================ + // Model constants specific to this subclass + // ================================================================ + + double McV = userInputs.get_model_constant_double("McV"); + double KcV = userInputs.get_model_constant_double("KcV"); + double WcV = userInputs.get_model_constant_double("WcV"); + double c0 = userInputs.get_model_constant_double("c0"); + double icamplitude = userInputs.get_model_constant_double("icamplitude"); + + // Declaring random number generator (Type std::mt19937_64) + engine rng; + // Declaring distribution (Type std::uniform_real_distribution) + distribution dist; + + // ================================================================ +}; diff --git a/automatic_tests/spinodalDecomposition/equations.cc b/automatic_tests/spinodalDecomposition/equations.cc new file mode 100644 index 000000000..467f47544 --- /dev/null +++ b/automatic_tests/spinodalDecomposition/equations.cc @@ -0,0 +1,130 @@ +// ================================================================================= +// Set the attributes of the primary field variables +// ================================================================================= +// This function sets attributes for each variable/equation in the app. The +// attributes are set via standardized function calls. The first parameter for +// each function call is the variable index (starting at zero). The first set of +// variable/equation attributes are the variable name (any string), the variable +// type (SCALAR/VECTOR), and the equation type (EXPLICIT_TIME_DEPENDENT/ +// TIME_INDEPENDENT/AUXILIARY). The next set of attributes describe the +// dependencies for the governing equation on the values and derivatives of the +// other variables for the value term and gradient term of the RHS and the LHS. +// The final pair of attributes determine whether a variable represents a field +// that can nucleate and whether the value of the field is needed for nucleation +// rate calculations. + +void +variableAttributeLoader::loadVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "c"); + set_variable_type(0, SCALAR); + set_variable_equation_type(0, EXPLICIT_TIME_DEPENDENT); + + set_dependencies_value_term_RHS(0, "c"); + set_dependencies_gradient_term_RHS(0, "grad(mu)"); + + // Variable 1 + set_variable_name(1, "mu"); + set_variable_type(1, SCALAR); + set_variable_equation_type(1, AUXILIARY); + + set_dependencies_value_term_RHS(1, "c"); + set_dependencies_gradient_term_RHS(1, "grad(c)"); +} + +// ============================================================================================= +// explicitEquationRHS (needed only if one or more equation is explict time +// dependent) +// ============================================================================================= +// This function calculates the right-hand-side of the explicit time-dependent +// equations for each variable. It takes "variable_list" as an input, which is a +// list of the value and derivatives of each of the variables at a specific +// quadrature point. The (x,y,z) location of that quadrature point is given by +// "q_point_loc". The function outputs two terms to variable_list -- one +// proportional to the test function and one proportional to the gradient of the +// test function. The index for each variable in this list corresponds to the +// index given at the top of this file. + +template +void +customPDE::explicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + scalarvalueType c = variable_list.get_scalar_value(0); + scalargradType mux = variable_list.get_scalar_gradient(1); + + // --- Setting the expressions for the terms in the governing equations --- + scalarvalueType eq_c = c; + scalargradType eqx_c = constV(-McV * userInputs.dtValue) * mux; + + // --- Submitting the terms for the governing equations --- + variable_list.set_scalar_value_term_RHS(0, eq_c); + variable_list.set_scalar_gradient_term_RHS(0, eqx_c); +} + +// ============================================================================================= +// nonExplicitEquationRHS (needed only if one or more equation is time +// independent or auxiliary) +// ============================================================================================= +// This function calculates the right-hand-side of all of the equations that are +// not explicit time-dependent equations. It takes "variable_list" as an input, +// which is a list of the value and derivatives of each of the variables at a +// specific quadrature point. The (x,y,z) location of that quadrature point is +// given by "q_point_loc". The function outputs two terms to variable_list -- +// one proportional to the test function and one proportional to the gradient of +// the test function. The index for each variable in this list corresponds to +// the index given at the top of this file. + +template +void +customPDE::nonExplicitEquationRHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + scalarvalueType c = variable_list.get_scalar_value(0); + scalargradType cx = variable_list.get_scalar_gradient(0); + + // --- Setting the expressions for the terms in the governing equations --- + + // The derivative of the local free energy + scalarvalueType fcV = 1.0 * WcV * c * (c - 1.0) * (c - 0.5); + + // The terms for the governing equations + scalarvalueType eq_mu = fcV; + scalargradType eqx_mu = constV(KcV) * cx; + + // --- Submitting the terms for the governing equations --- + + variable_list.set_scalar_value_term_RHS(1, eq_mu); + variable_list.set_scalar_gradient_term_RHS(1, eqx_mu); +} + +// ============================================================================================= +// equationLHS (needed only if at least one equation is time independent) +// ============================================================================================= +// This function calculates the left-hand-side of time-independent equations. It +// takes "variable_list" as an input, which is a list of the value and +// derivatives of each of the variables at a specific quadrature point. The +// (x,y,z) location of that quadrature point is given by "q_point_loc". The +// function outputs two terms to variable_list -- one proportional to the test +// function and one proportional to the gradient of the test function -- for the +// left-hand-side of the equation. The index for each variable in this list +// corresponds to the index given at the top of this file. If there are multiple +// elliptic equations, conditional statements should be sed to ensure that the +// correct residual is being submitted. The index of the field being solved can +// be accessed by "this->currentFieldIndex". + +template +void +customPDE::equationLHS( + [[maybe_unused]] variableContainer> &variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{} diff --git a/automatic_tests/spinodalDecomposition/gold_integratedFields.txt b/automatic_tests/spinodalDecomposition/gold_integratedFields.txt new file mode 100644 index 000000000..5845b63e7 --- /dev/null +++ b/automatic_tests/spinodalDecomposition/gold_integratedFields.txt @@ -0,0 +1,11 @@ +0 f_tot 40.9682267 +4 f_tot 37.66188071 +8 f_tot 30.9047905 +12 f_tot 27.4978036 +16 f_tot 24.88328855 +20 f_tot 22.87153383 +24 f_tot 21.45111207 +28 f_tot 20.38772127 +32 f_tot 19.53766885 +36 f_tot 18.82222235 +40 f_tot 18.25004626 diff --git a/automatic_tests/spinodalDecomposition/parameters.prm b/automatic_tests/spinodalDecomposition/parameters.prm new file mode 100644 index 000000000..96c06c107 --- /dev/null +++ b/automatic_tests/spinodalDecomposition/parameters.prm @@ -0,0 +1,80 @@ +# ================================================================================= +# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) +# ================================================================================= +set Number of dimensions = 2 + +# ================================================================================= +# Set the length of the domain in all three dimensions +# (Domain size Z ignored in 2D) +# ================================================================================= +set Domain size X = 51.2 +set Domain size Y = 51.2 +set Domain size Z = 100 + +# ================================================================================= +# Set the element parameters +# ================================================================================= +set Subdivisions X = 1 +set Subdivisions Y = 1 +set Subdivisions Z = 1 + +set Refine factor = 8 + +set Element degree = 1 + +# ================================================================================= +# Set the adaptive mesh refinement parameters +# ================================================================================= +set Mesh adaptivity = true + +set Max refinement level = 8 +set Min refinement level = 2 + +set Steps between remeshing operations = 2000 + +subsection Refinement criterion: c + set Criterion type = VALUE + set Value lower bound = 0.01 + set Value upper bound = 0.99 +end + +# ================================================================================= +# Set the time step parameters +# ================================================================================= +set Time step = 0.004 + +set Number of time steps = 10000 + +# ================================================================================= +# Set the boundary conditions +# ================================================================================= +set Boundary condition for variable c = PERIODIC +set Boundary condition for variable mu = PERIODIC + +# ================================================================================= +# Set the model constants +# ================================================================================= + +# The mobility, McV in equations.h +set Model constant McV = 1.0, DOUBLE + +# The gradient energy coefficient, KcV in equations.h +set Model constant KcV = 0.01, DOUBLE + +# The double well coefficient +set Model constant WcV = 1.0, DOUBLE + +# The average composition +set Model constant c0 = 0.50, DOUBLE + +# The initial perturbation amplitude +set Model constant icamplitude = 0.01, DOUBLE + +# ================================================================================= +# Set the output parameters +# ================================================================================= +set Output condition = EQUAL_SPACING + +set Number of outputs = 10 + +set Skip print steps = 1000 diff --git a/automatic_tests/spinodalDecomposition/postprocess.cc b/automatic_tests/spinodalDecomposition/postprocess.cc new file mode 100644 index 000000000..d4486ebc6 --- /dev/null +++ b/automatic_tests/spinodalDecomposition/postprocess.cc @@ -0,0 +1,78 @@ +// ============================================================================================= +// loadPostProcessorVariableAttributes: Set the attributes of the postprocessing +// variables +// ============================================================================================= +// This function is analogous to 'loadVariableAttributes' in 'equations.h', but +// for the postprocessing expressions. It sets the attributes for each +// postprocessing expression, including its name, whether it is a vector or +// scalar (only scalars are supported at present), its dependencies on other +// variables and their derivatives, and whether to calculate an integral of the +// postprocessed quantity over the entire domain. Note: this function is not a +// member of customPDE. + +void +variableAttributeLoader::loadPostProcessorVariableAttributes() +{ + // Variable 0 + set_variable_name(0, "f_tot"); + set_variable_type(0, SCALAR); + + set_dependencies_value_term_RHS(0, "c, grad(c)"); + set_dependencies_gradient_term_RHS(0, ""); + + set_output_integral(0, true); +} + +// ============================================================================================= +// postProcessedFields: Set the postprocessing expressions +// ============================================================================================= +// This function is analogous to 'explicitEquationRHS' and +// 'nonExplicitEquationRHS' in equations.h. It takes in "variable_list" and +// "q_point_loc" as inputs and outputs two terms in the expression for the +// postprocessing variable -- one proportional to the test function and one +// proportional to the gradient of the test function. The index for each +// variable in this list corresponds to the index given at the top of this file +// (for submitting the terms) and the index in 'equations.h' for assigning the +// values/derivatives of the primary variables. + +template +void +customPDE::postProcessedFields( + [[maybe_unused]] const variableContainer> + &variable_list, + [[maybe_unused]] variableContainer> + &pp_variable_list, + [[maybe_unused]] const Point> q_point_loc, + [[maybe_unused]] const VectorizedArray element_volume) const +{ + // --- Getting the values and derivatives of the model variables --- + + // The concentration and its derivatives + scalarvalueType c = variable_list.get_scalar_value(0); + scalargradType cx = variable_list.get_scalar_gradient(0); + + // --- Setting the expressions for the terms in the postprocessing expressions + // --- + + scalarvalueType f_tot = constV(0.0); + + // The homogenous free energy + scalarvalueType f_chem = 0.25 * WcV * (c * c * c * c - 2.0 * c * c * c + c * c); + + // The gradient free energy + scalarvalueType f_grad = constV(0.0); + + for (int i = 0; i < dim; i++) + { + for (int j = 0; j < dim; j++) + { + f_grad += constV(0.5 * KcV) * cx[i] * cx[j]; + } + } + + // The total free energy + f_tot = f_chem + f_grad; + + // --- Submitting the terms for the postprocessing expressions --- + pp_variable_list.set_scalar_value_term_RHS(0, f_tot); +} diff --git a/tests/automatic_tests/test_results.txt b/automatic_tests/test_results.txt similarity index 97% rename from tests/automatic_tests/test_results.txt rename to automatic_tests/test_results.txt index 31a30082b..455d572dd 100644 --- a/tests/automatic_tests/test_results.txt +++ b/automatic_tests/test_results.txt @@ -1792,3 +1792,54 @@ Time: 66.96555471420288 Tests Passed: 5/5 --------------------------------------------------------- +--------------------------------------------------------- +Regression test on 2024-12-11 12:27 +Architecture: x86_64 +Model name: AMD Ryzen 9 7900X 12-Core Processor +CPU(s): 24 +CPU max/min MHz: 4701.0000, 0.0000 +Hypervisor vendor: Windows Subsystem for Linux +Number of processes: 1 +--------------------------------------------------------- +Application: allenCahn +Result: Pass +Time: 4.505462408065796 + +Application: cahnHilliard +Result: Pass +Time: 2.7520039081573486 + +Application: CHAC_anisotropyRegularized +Result: Pass +Time: 4.074458122253418 + +Application: CHiMaD_benchmark6a +Result: New Gold Standard +Time: 23.678637266159058 + +Application: corrosion_microgalvanic +Result: New Gold Standard +Time: 20.824553728103638 + +Application: coupledCahnHilliardAllenCahn +Result: Pass +Time: 2.059581995010376 + +Application: grainGrowth +Result: New Gold Standard +Time: 12.296957015991211 + +Application: precipitateEvolution +Result: Pass +Time: 14.774498701095581 + +Application: precipitateEvolution_pfunction +Result: New Gold Standard +Time: 21.35426139831543 + +Application: spinodalDecomposition +Result: New Gold Standard +Time: 23.352713108062744 + +Tests Passed: 10/10 +--------------------------------------------------------- diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 000000000..8aec4ac4d --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,125 @@ +## +# CMake for the PRISMS-PF unit tests +# Adapted from the ASPECT CMake file +## + +cmake_minimum_required(VERSION 3.1.0) + +project(prisms_pf_unit_tests) + +# ========================================================= +# Some basic bookkeeping +# ========================================================= + +# Check that a prior CMakeCache is not located in the build directory +if(EXISTS ${CMAKE_SOURCE_DIR}/CMakeCache.txt) + message(FATAL_ERROR "Detected the file:\n" + "${CMAKE_SOURCE_DIR}/CMakeCache.txt\n" + "in your source directory, which may be leftover from prior builds. " + "Please delete the file before running cmake again.") +endif() + +# ========================================================= +# External libraries +# ========================================================= + +# Find deal.II installation +find_package(deal.II 9.2.0 QUIET + HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} + ) +if(NOT ${deal.II_FOUND}) + message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " + "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " + "recent version of deal.II." + ) +endif() + +message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") + +set(DEALII_INSTALL_VALID ON) + +if(NOT DEAL_II_WITH_P4EST) + message(SEND_ERROR + "\n**deal.II was built without support for p4est!\n" + ) + set(DEALII_INSTALL_VALID OFF) +endif() + +if(NOT DEALII_INSTALL_VALID) + message(FATAL_ERROR + "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" + ) +endif() + +deal_ii_initialize_cached_variables() + +# Create compile_commands.json +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +set(FORCE_COLORED_OUTPUT ON CACHE BOOL "Forces colored output when compiling with gcc and clang.") + +# Setting up tests +MESSAGE(STATUS "Setting up tests with CMAKE_BUILD_TYPE ${CMAKE_BUILD_TYPE}") + +enable_testing() + +add_custom_target(tests) + +# diff tool +find_program(DIFF_EXECUTABLE + NAMES diff + HINTS ${DIFF_DIR} + PATH_SUFFIXES bin + ) +find_program(NUMDIFF_EXECUTABLE + NAMES numdiff + HINTS ${NUMDIFF_DIR} + PATH_SUFFIXES bin + ) + +mark_as_advanced(DIFF_EXECUTABLE NUMDIFF_EXECUTABLE) + +if("${TEST_DIFF}" STREQUAL "") + if(NOT NUMDIFF_EXECUTABLE MATCHES "-NOTFOUND") + set(TEST_DIFF ${NUMDIFF_EXECUTABLE}) + elseif(NOT DIFF_EXECUTABLE MATCHES "-NOTFOUND") + set(TEST_DIFF ${DIFF_EXECUTABLE}) + else() + message(FATAL_ERROR + "Could not find diff or numdiff. One of those must be installed for running the testsuite./n" + "Please specify TEST_DIFF by hand." + ) + endif() +endif() + +# Set the name of the project and target: +set(TARGET "main") + +file(GLOB_RECURSE TEST_SOURCES RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} "*.cc") + +# Set location of Catch +include_directories(${CMAKE_SOURCE_DIR}/../contrib/catch/) + +# Set location of core library files +include_directories(${CMAKE_SOURCE_DIR}/../include) + +# Declare all source files the target consists of: +set(TARGET_SRC + main.cc + ${TEST_SOURCES} + ) + +# Check if there has been updates to main library +set(dir ${PROJECT_SOURCE_DIR}/../..) +EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) +EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) + +# Set targets & link libraries for the build type +add_executable(main ${TARGET_SRC}) +set_property(TARGET main PROPERTY OUTPUT_NAME main) +deal_ii_setup_target(main DEBUG) +target_link_libraries(main ${CMAKE_SOURCE_DIR}/../libprisms-pf-debug.a) + +# CTest +add_test(NAME PRISMS_PF_Testsuite COMMAND main) \ No newline at end of file diff --git a/tests/automatic_tests/CHAC_anisotropyRegularized/parameters.prm b/tests/automatic_tests/CHAC_anisotropyRegularized/parameters.prm deleted file mode 100644 index dcf513c24..000000000 --- a/tests/automatic_tests/CHAC_anisotropyRegularized/parameters.prm +++ /dev/null @@ -1,123 +0,0 @@ -# ================================================================================= -# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) -# ================================================================================= -set Number of dimensions = 2 - -# ================================================================================= -# Set the length of the domain in all three dimensions -# (Domain size Z ignored in 2D) -# ================================================================================= -# Each axes spans from zero to the specified length -set Domain size X = 100 -set Domain size Y = 100 -set Domain size Z = 100 - -# ================================================================================= -# Set the element parameters -# ================================================================================= -# The number of elements in each direction is 2^(refineFactor) * subdivisions -# Subdivisions Z ignored in 2D -# For optimal performance, use refineFactor primarily to determine the element size -set Subdivisions X = 1 -set Subdivisions Y = 1 -set Subdivisions Z = 1 - -set Refine factor = 7 - -# Set the polynomial degree of the element (allowed values: 1, 2, or 3) -set Element degree = 1 - -# ================================================================================= -# Set the adaptive mesh refinement parameters -# ================================================================================= -# Set the flag determining if adaptive meshing is activated -set Mesh adaptivity = true - -# Set the maximum and minimum level of refinement -# When adaptive meshing is enabled, the refine factor set in the block above is -# only used to generate the first pass of the mesh as the initial conditions are -# applied. It should be set somewhere between the max and min levels below. -set Max refinement level = 7 -set Min refinement level = 4 - -# Set the number of time steps between remeshing operations -set Steps between remeshing operations = 1000 - -# Set the criteria for adapting the mesh -subsection Refinement criterion: n - # Select whether the mesh is refined based on the variable value (VALUE), - # its gradient (GRADIENT), or both (VALUE_AND_GRADIENT) - set Criterion type = VALUE - # Set the lower and upper bounds for the value-based refinement window - set Value lower bound = 0.001 - set Value upper bound = 0.999 -end - -# ================================================================================= -# Set the time step parameters -# ================================================================================= -# The size of the time step -set Time step = 5.0e-2 - -# The simulation ends when either the number of time steps is reached or the -# simulation time is reached. -set Number of time steps = 5000 - -# ================================================================================= -# Set the boundary conditions -# ================================================================================= -# Set the boundary condition for each variable, where each variable is given by -# its name, as defined in equations.h. The four boundary condition -# types are NATURAL, DIRICHLET, NON_UNIFORM_DIRICHLET and PERIODIC. If all -# of the boundaries have the same boundary condition, only one boundary condition -# type needs to be given. If multiple boundary condition types are needed, give a -# comma-separated list of the types. The order is the miniumum of x, maximum of x, -# minimum of y, maximum of y, minimum of z, maximum of z (i.e left, right, bottom, -# top in 2D and left, right, bottom, top, front, back in 3D). The value of a -# Dirichlet BC is specfied in the following way -- DIRCHILET: val -- where 'val' is -# the desired value. If the boundary condition is NON_UNIFORM_DIRICHLET, the -# boundary condition should be specified in the appropriate function in 'ICs_and_BCs.h'. -# Example 1: All periodic BCs for variable 'c' -# set Boundary condition for variable c = PERIODIC -# Example 2: Zero-derivative BCs on the left and right, Dirichlet BCs with value -# 1.5 on the top and bottom for variable 'n' in 2D -# set Boundary condition for variable n = NATURAL, NATURAL, DIRICHLET: 1.5, DIRICHLET: 1.5 - -set Boundary condition for variable c = NATURAL -set Boundary condition for variable n = NATURAL -set Boundary condition for variable biharm = NATURAL - -# ================================================================================= -# Set the model constants -# ================================================================================= -# Set the user-defined model constants, which must have a counter-part given in -# customPDE.h. These are most often used in the residual equations in equations.h, -# but may also be used for initial conditions and nucleation calculations. The type -# options currently are DOUBLE, INT, BOOL, TENSOR, and [symmetry] ELASTIC CONSTANTS -# where [symmetry] is ISOTROPIC, TRANSVERSE, ORTHOTROPIC, or ANISOTROPIC. - -# The CH mobility, McV in equations.h -set Model constant McV = 1.0, DOUBLE - -# The AC mobility, MnV in equations.h -set Model constant MnV = 0.1, DOUBLE - -# Anisotropy parameter -set Model constant epsilonM = 0.2, DOUBLE - -# Regularization parameter -set Model constant delta2 = 1.0, DOUBLE - -# ================================================================================= -# Set the output parameters -# ================================================================================= -# Type of spacing between outputs ("EQUAL_SPACING", "LOG_SPACING", "N_PER_DECADE", -# or "LIST") -set Output condition = EQUAL_SPACING - -# Number of times the program outputs the fields (total number for "EQUAL_SPACING" -# and "LOG_SPACING", number per decade for "N_PER_DECADE", ignored for "LIST") -set Number of outputs = 10 - -# The number of time steps between updates being printed to the screen -set Skip print steps = 1000 diff --git a/tests/automatic_tests/allenCahn/parameters.prm b/tests/automatic_tests/allenCahn/parameters.prm deleted file mode 100644 index 6fa3f6433..000000000 --- a/tests/automatic_tests/allenCahn/parameters.prm +++ /dev/null @@ -1,92 +0,0 @@ -# ================================================================================= -# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) -# ================================================================================= -set Number of dimensions = 2 - -# ================================================================================= -# Set the length of the domain in all three dimensions -# (Domain size Z ignored in 2D) -# ================================================================================= -# Each axes spans from zero to the specified length -set Domain size X = 100 -set Domain size Y = 100 -set Domain size Z = 100 - -# ================================================================================= -# Set the element parameters -# ================================================================================= -# The number of elements in each direction is 2^(refineFactor) * subdivisions -# Subdivisions Z ignored in 2D -# For optimal performance, use refineFactor primarily to determine the element size -set Subdivisions X = 1 -set Subdivisions Y = 1 -set Subdivisions Z = 1 - -set Refine factor = 8 - -# Set the polynomial degree of the element (allowed values: 1, 2, or 3) -set Element degree = 1 - -# ================================================================================= -# Set the time step parameters -# ================================================================================= -# The size of the time step -set Time step = 1.0e-2 - -# The simulation ends when either the number of time steps is reached or the -# simulation time is reached. -set Number of time steps = 5000 - -# ================================================================================= -# Set the output parameters -# ================================================================================= -# Type of spacing between outputs ("EQUAL_SPACING", "LOG_SPACING", "N_PER_DECADE", -# or "LIST") -set Output condition = EQUAL_SPACING - -# Number of times the program outputs the fields (total number for "EQUAL_SPACING" -# and "LOG_SPACING", number per decade for "N_PER_DECADE", ignored for "LIST") -set Number of outputs = 5 - -# Whether to print timing information every time the code outputs -set Print timing information with output = true - -# The number of time steps between updates being printed to the screen -set Skip print steps = 1000 - -# ================================================================================= -# Set the boundary conditions -# ================================================================================= -# Set the boundary condition for each variable, where each variable is given by -# its name, as defined in equations.cc. The four boundary condition -# types are NATURAL, DIRICHLET, NON_UNIFORM_DIRICHLET and PERIODIC. If all -# of the boundaries have the same boundary condition, only one boundary condition -# type needs to be given. If multiple boundary condition types are needed, give a -# comma-separated list of the types. The order is the miniumum of x, maximum of x, -# minimum of y, maximum of y, minimum of z, maximum of z (i.e left, right, bottom, -# top in 2D and left, right, bottom, top, front, back in 3D). The value of a -# Dirichlet BC is specfied in the following way -- DIRCHILET: val -- where 'val' is -# the desired value. If the boundary condition is NON_UNIFORM_DIRICHLET, the -# boundary condition should be specified in the appropriate function in 'ICs_and_BCs.h'. -# Example 1: All periodic BCs for variable 'c' -# set Boundary condition for variable c = PERIODIC -# Example 2: Zero-derivative BCs on the left and right, Dirichlet BCs with value -# 1.5 on the top and bottom for variable 'n' in 2D -# set Boundary condition for variable n = NATURAL, NATURAL, DIRICHLET: 1.5, DIRICHLET: 1.5 - -set Boundary condition for variable n = NATURAL - -# ================================================================================= -# Set the model constants -# ================================================================================= -# Set the user-defined model constants, which must have a counter-part given in -# customPDE.h. These are most often used in the residual equations in equations.cc, -# but may also be used for initial conditions and nucleation calculations. The type -# options currently are DOUBLE, INT, BOOL, TENSOR, and [symmetry] ELASTIC CONSTANTS -# where [symmetry] is ISOTROPIC, TRANSVERSE, ORTHOTROPIC, or ANISOTROPIC. - -# The mobility, MnV in equations.cc -set Model constant MnV = 1.0, DOUBLE - -# The gradient energy coefficient, KnV in equations.cc -set Model constant KnV = 2.0, DOUBLE diff --git a/tests/automatic_tests/cahnHilliard/parameters.prm b/tests/automatic_tests/cahnHilliard/parameters.prm deleted file mode 100644 index a463d5567..000000000 --- a/tests/automatic_tests/cahnHilliard/parameters.prm +++ /dev/null @@ -1,90 +0,0 @@ -# ================================================================================= -# Set the number of dimensions (2 or 3 for a 2D or 3D calculation) -# ================================================================================= -set Number of dimensions = 2 - -# ================================================================================= -# Set the length of the domain in all three dimensions -# (Domain size Z ignored in 2D) -# ================================================================================= -# Each axes spans from zero to the specified length -set Domain size X = 100 -set Domain size Y = 100 -set Domain size Z = 100 - -# ================================================================================= -# Set the element parameters -# ================================================================================= -# The number of elements in each direction is 2^(refineFactor) * subdivisions -# Subdivisions Z ignored in 2D -# For optimal performance, use refineFactor primarily to determine the element size -set Subdivisions X = 1 -set Subdivisions Y = 1 -set Subdivisions Z = 1 - -set Refine factor = 6 - -# Set the polynomial degree of the element (allowed values: 1, 2, or 3) -set Element degree = 2 - -# ================================================================================= -# Set the time step parameters -# ================================================================================= -# The size of the time step -set Time step = 1.0e-3 - -# The simulation ends when either the number of time steps is reached or the -# simulation time is reached. -set Number of time steps = 5000 - -# ================================================================================= -# Set the boundary conditions -# ================================================================================= -# Set the boundary condition for each variable, where each variable is given by -# its name, as defined in equations.h. The four boundary condition -# types are NATURAL, DIRICHLET, NON_UNIFORM_DIRICHLET and PERIODIC. If all -# of the boundaries have the same boundary condition, only one boundary condition -# type needs to be given. If multiple boundary condition types are needed, give a -# comma-separated list of the types. The order is the miniumum of x, maximum of x, -# minimum of y, maximum of y, minimum of z, maximum of z (i.e left, right, bottom, -# top in 2D and left, right, bottom, top, front, back in 3D). The value of a -# Dirichlet BC is specfied in the following way -- DIRCHILET: val -- where 'val' is -# the desired value. If the boundary condition is NON_UNIFORM_DIRICHLET, the -# boundary condition should be specified in the appropriate function in 'ICs_and_BCs.h'. -# Example 1: All periodic BCs for variable 'c' -# set Boundary condition for variable c = PERIODIC -# Example 2: Zero-derivative BCs on the left and right, Dirichlet BCs with value -# 1.5 on the top and bottom for variable 'n' in 2D -# set Boundary condition for variable n = NATURAL, NATURAL, DIRICHLET: 1.5, DIRICHLET: 1.5 - -set Boundary condition for variable c = NATURAL -set Boundary condition for variable mu = NATURAL - -# ================================================================================= -# Set the model constants -# ================================================================================= -# Set the user-defined model constants, which must have a counter-part given in -# customPDE.h. These are most often used in the residual equations in equations.h, -# but may also be used for initial conditions and nucleation calculations. The type -# options currently are DOUBLE, INT, BOOL, TENSOR, and [symmetry] ELASTIC CONSTANTS -# where [symmetry] is ISOTROPIC, TRANSVERSE, ORTHOTROPIC, or ANISOTROPIC. - -# The mobility, McV in equations.h -set Model constant McV = 1.0, DOUBLE - -# The gradient energy coefficient, KcV in equations.h -set Model constant KcV = 1.5, DOUBLE - -# ================================================================================= -# Set the output parameters -# ================================================================================= -# Type of spacing between outputs ("EQUAL_SPACING", "LOG_SPACING", "N_PER_DECADE", -# or "LIST") -set Output condition = EQUAL_SPACING - -# Number of times the program outputs the fields (total number for "EQUAL_SPACING" -# and "LOG_SPACING", number per decade for "N_PER_DECADE", ignored for "LIST") -set Number of outputs = 10 - -# The number of time steps between updates being printed to the screen -set Skip print steps = 1000 diff --git a/tests/automatic_tests/precipitateEvolution/CMakeLists.txt b/tests/automatic_tests/precipitateEvolution/CMakeLists.txt deleted file mode 100644 index 3b7e5c118..000000000 --- a/tests/automatic_tests/precipitateEvolution/CMakeLists.txt +++ /dev/null @@ -1,165 +0,0 @@ -## -# CMake script for the PRISMS-PF applications -# Adapted from the ASPECT CMake file -## - -cmake_minimum_required(VERSION 3.1.0) - -project(myapp) - -message(STATUS "") -message(STATUS "=========================================================") -message(STATUS "Configuring PRISMS-PF application") -message(STATUS "=========================================================") -message(STATUS "") - -# Setup the build type (debug, release, debugrelease) -if("${CMAKE_BUILD_TYPE}" STREQUAL "") - set(CMAKE_BUILD_TYPE "DebugRelease" - CACHE STRING - "Choose the type of build, options are: Debug, Release and DebugRelease." - FORCE) -endif() - -# Convert build type into the debug and release builds, which may or may -# not be built. -if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR - "${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR - "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease" ) - message(STATUS "Setting up PRISMS-PF application for ${CMAKE_BUILD_TYPE} mode.") -else() - message(FATAL_ERROR - "CMAKE_BUILD_TYPE must either be 'Release', 'Debug', or 'DebugRelease', but is set to '${CMAKE_BUILD_TYPE}'.") -endif() - -if("${CMAKE_BUILD_TYPE}" STREQUAL "Debug" OR - "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") - set(PRISMS_PF_BUILD_DEBUG "ON") -else() - set(PRISMS_PF_BUILD_DEBUG "OFF") -endif() - -if("${CMAKE_BUILD_TYPE}" STREQUAL "Release" OR - "${CMAKE_BUILD_TYPE}" STREQUAL "DebugRelease") - set(PRISMS_PF_BUILD_RELEASE "ON") -else() - set(PRISMS_PF_BUILD_RELEASE "OFF") -endif() - -# Find deal.II installation -find_package(deal.II 9.2.0 QUIET - HINTS ${DEAL_II_DIR} $ENV{DEAL_II_DIR} - ) -if(NOT ${deal.II_FOUND}) - message(FATAL_ERROR "\n*** Could not find a recent version of deal.II. ***\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake " - "or set an environment variable \"DEAL_II_DIR\" that contains a path to a " - "recent version of deal.II." - ) -endif() - -message(STATUS "Found deal.II version ${DEAL_II_PACKAGE_VERSION} at '${deal.II_DIR}'") - -set(DEALII_INSTALL_VALID ON) - -if(NOT DEAL_II_WITH_P4EST) - message(SEND_ERROR - "\n**deal.II was built without support for p4est!\n" - ) - set(DEALII_INSTALL_VALID OFF) -endif() - -if(NOT DEALII_INSTALL_VALID) - message(FATAL_ERROR - "\nPRISMS-PF requires a deal.II installation with certain features enabled!\n" - ) -endif() - -deal_ii_initialize_cached_variables() - -# Make and ninja build options -if(CMAKE_GENERATOR MATCHES "Ninja") - set(_make_command "$ ninja") -else() - set(_make_command " $ make") -endif() - -# Debug and release targets -if(${DEAL_II_BUILD_TYPE} MATCHES "DebugRelease") - add_custom_target(release - COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Release . - COMMAND ${CMAKE_COMMAND} -E echo "***" - COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Release mode. Now recompile with: ${_make_command}" - COMMAND ${CMAKE_COMMAND} -E echo "***" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - VERBATIM - COMMENT "switching to RELEASE mode..." - ) - - add_custom_target(debug - COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=Debug . - COMMAND ${CMAKE_COMMAND} -E echo "***" - COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug mode. Now recompile with: ${_make_command}" - COMMAND ${CMAKE_COMMAND} -E echo "***" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - VERBATIM - COMMENT "switching to DEBUG mode..." - ) - - add_custom_target(debugrelease - COMMAND ${CMAKE_COMMAND} -D CMAKE_BUILD_TYPE=DebugRelease . - COMMAND ${CMAKE_COMMAND} -E echo "***" - COMMAND ${CMAKE_COMMAND} -E echo "*** Switched to Debug and Release mode. Now recompile with: ${_make_command}" - COMMAND ${CMAKE_COMMAND} -E echo "***" - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - VERBATIM - COMMENT "switching to DEBUG/RELEASE mode..." - ) -endif() - -# Add distclean target to clean build -add_custom_target(distclean - COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target clean - COMMAND ${CMAKE_COMMAND} -E remove_directory CMakeFiles - COMMAND ${CMAKE_COMMAND} -E remove - CMakeCache.txt cmake_install.cmake Makefile - build.ninja rules.ninja .ninja_deps .ninja_log - COMMENT "distclean invoked" - ) - -# Add postprocess.cc and nucleation.cc if they exist -if(EXISTS "postprocess.cc") - add_definitions(-DPOSTPROCESS_FILE_EXISTS) -endif() -if(EXISTS "nucleation.cc") - add_definitions(-DNUCLEATION_FILE_EXISTS) -endif() - -# Set location of files -include_directories(${CMAKE_SOURCE_DIR}/../../../include) -include_directories(${CMAKE_SOURCE_DIR}/../../../src) -include_directories(${CMAKE_SOURCE_DIR}) - -# Set the location of the main.cc file -set(TARGET_SRC "${CMAKE_SOURCE_DIR}/../main.cc") - -# Check if there has been updates to main library -set(dir ${PROJECT_SOURCE_DIR}/../../..) -EXECUTE_PROCESS(COMMAND "rm" "CMakeCache.txt" WORKING_DIRECTORY ${dir}) -EXECUTE_PROCESS(COMMAND "cmake" "CMakeLists.txt" WORKING_DIRECTORY ${dir}) -EXECUTE_PROCESS(COMMAND "make" WORKING_DIRECTORY ${dir}) - -# Set targets & link libraries for the build type -if(${PRISMS_PF_BUILD_DEBUG} STREQUAL "ON") - add_executable(main_debug ${TARGET_SRC}) - set_property(TARGET main_debug PROPERTY OUTPUT_NAME main-debug) - deal_ii_setup_target(main_debug DEBUG) - target_link_libraries(main_debug ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-debug.a) -endif() - -if(${PRISMS_PF_BUILD_RELEASE} STREQUAL "ON") - add_executable(main_release ${TARGET_SRC}) - set_property(TARGET main_release PROPERTY OUTPUT_NAME main) - deal_ii_setup_target(main_release RELEASE) - target_link_libraries(main_release ${CMAKE_SOURCE_DIR}/../../../libprisms-pf-release.a) -endif() \ No newline at end of file diff --git a/tests/core/boundary_conditions/boundary_conditions.cc b/tests/core/boundary_conditions/boundary_conditions.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/core/compute_invM.cc b/tests/core/compute_invM.cc new file mode 100644 index 000000000..2e8b76cb3 --- /dev/null +++ b/tests/core/compute_invM.cc @@ -0,0 +1,47 @@ +#include +#include +#include +#include +#include +#include + +#include "catch.hpp" + +// #include +#include +#include + +TEST_CASE("Compute invM") +{ + SECTION("1D Subdivided Hyper Rectangle") + { + SECTION("Scalar 1st Degree") + { + unsigned int dim = 1; + unsigned int degree = 1; + + // Initialize triangulation for a unit cube + std::vector subdivisions(dim, 1); + + dealii::Triangulation<1> triangulation; + dealii::GridGenerator::subdivided_hyper_rectangle(triangulation, + subdivisions, + dealii::Point<1>(), + dealii::Point<1>(1.0)); + + // Create scalar finite element space + dealii::FE_Q<1> fe_q(dealii::QGaussLobatto<1>(degree + 1)); + dealii::FESystem<1> fe_system(fe_q, 1); + + // Create DoFHandler + dealii::DoFHandler<1> dof_handler(triangulation); + dof_handler.distribute_dofs(fe_system); + + REQUIRE(dof_handler.n_dofs() == 2); + } + } + SECTION("2D") + {} + SECTION("3D") + {} +} \ No newline at end of file diff --git a/tests/core/initial_conditions/initial_conditions.cc b/tests/core/initial_conditions/initial_conditions.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/core/postprocessing/postprocessing.cc b/tests/core/postprocessing/postprocessing.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/core/refinement/refinement.cc b/tests/core/refinement/refinement.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/core/solvers/solvers.cc b/tests/core/solvers/solvers.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/field_input/field_input.cc b/tests/field_input/field_input.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/grains/grains.cc b/tests/grains/grains.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/main.cc b/tests/main.cc new file mode 100644 index 000000000..063e87874 --- /dev/null +++ b/tests/main.cc @@ -0,0 +1,2 @@ +#define CATCH_CONFIG_MAIN +#include "catch.hpp" \ No newline at end of file diff --git a/tests/nucleation/nucleation.cc b/tests/nucleation/nucleation.cc new file mode 100644 index 000000000..e69de29bb diff --git a/tests/utilities/utilities.cc b/tests/utilities/utilities.cc new file mode 100644 index 000000000..e69de29bb