From a1da5e8b1c5db570b26d7a9e819dab04d58788c6 Mon Sep 17 00:00:00 2001 From: aradermacher Date: Mon, 18 Sep 2023 11:27:13 +0200 Subject: [PATCH] add parameter dimensionalty test in material problems --- .../experimental_setup/base_experiment.py | 1 + .../experimental_setup/cantilever_beam.py | 2 +- .../experimental_setup/simple_beam.py | 1 - .../experimental_setup/simple_cube.py | 2 +- .../experimental_setup/tensile_beam.py | 2 +- .../finite_element_problem/base_material.py | 48 +++++++++++++++---- .../finite_element_problem/concrete_am.py | 22 +++++---- .../concrete_thermo_mechanical.py | 21 ++++---- .../linear_elasticity.py | 42 +++++++++++++--- .../test_linear_simple_cube.py | 5 ++ .../test_thermo_mechanical_cube.py | 2 +- 11 files changed, 110 insertions(+), 38 deletions(-) diff --git a/src/fenicsxconcrete/experimental_setup/base_experiment.py b/src/fenicsxconcrete/experimental_setup/base_experiment.py index fe95ba5..c2974ee 100644 --- a/src/fenicsxconcrete/experimental_setup/base_experiment.py +++ b/src/fenicsxconcrete/experimental_setup/base_experiment.py @@ -39,6 +39,7 @@ def __init__(self, parameters: dict[str, pint.Quantity]) -> None: setup_parameters.update(parameters) # get logger info which parameters are set to default values + # plus check dimensionality of input parameters keys_set_default = [] for key in dict(default_p): if key not in parameters: diff --git a/src/fenicsxconcrete/experimental_setup/cantilever_beam.py b/src/fenicsxconcrete/experimental_setup/cantilever_beam.py index 7d2607e..7b06b38 100644 --- a/src/fenicsxconcrete/experimental_setup/cantilever_beam.py +++ b/src/fenicsxconcrete/experimental_setup/cantilever_beam.py @@ -70,7 +70,7 @@ def default_parameters() -> dict[str, pint.Quantity]: """ setup_parameters = {} - setup_parameters["degree"] = 2 * ureg("") # polynomial degree + setup_parameters["length"] = 1 * ureg("m") setup_parameters["height"] = 0.3 * ureg("m") setup_parameters["width"] = 0.3 * ureg("m") # only relevant for 3D case diff --git a/src/fenicsxconcrete/experimental_setup/simple_beam.py b/src/fenicsxconcrete/experimental_setup/simple_beam.py index ab2694f..edfecdb 100644 --- a/src/fenicsxconcrete/experimental_setup/simple_beam.py +++ b/src/fenicsxconcrete/experimental_setup/simple_beam.py @@ -74,7 +74,6 @@ def default_parameters() -> dict[str, pint.Quantity]: """ setup_parameters = {} - setup_parameters["degree"] = 2 * ureg("") # polynomial degree setup_parameters["load"] = 10000 * ureg("N/m^2") setup_parameters["length"] = 1 * ureg("m") diff --git a/src/fenicsxconcrete/experimental_setup/simple_cube.py b/src/fenicsxconcrete/experimental_setup/simple_cube.py index 463dcad..33e0fca 100644 --- a/src/fenicsxconcrete/experimental_setup/simple_cube.py +++ b/src/fenicsxconcrete/experimental_setup/simple_cube.py @@ -45,7 +45,7 @@ def default_parameters() -> dict[str, pint.Quantity]: """ setup_parameters = {} - setup_parameters["degree"] = 2 * ureg("") # polynomial degree + setup_parameters["height"] = 1 * ureg("m") setup_parameters["width"] = 1 * ureg("m") setup_parameters["length"] = 1 * ureg("m") diff --git a/src/fenicsxconcrete/experimental_setup/tensile_beam.py b/src/fenicsxconcrete/experimental_setup/tensile_beam.py index eefc92e..1567ef3 100644 --- a/src/fenicsxconcrete/experimental_setup/tensile_beam.py +++ b/src/fenicsxconcrete/experimental_setup/tensile_beam.py @@ -70,7 +70,7 @@ def default_parameters() -> dict[str, pint.Quantity]: """ setup_parameters = {} - setup_parameters["degree"] = 2 * ureg("") # polynomial degree + setup_parameters["length"] = 1 * ureg("m") setup_parameters["height"] = 0.3 * ureg("m") setup_parameters["width"] = 0.3 * ureg("m") # only relevant for 3D case diff --git a/src/fenicsxconcrete/finite_element_problem/base_material.py b/src/fenicsxconcrete/finite_element_problem/base_material.py index 6411b40..9d3f39e 100644 --- a/src/fenicsxconcrete/finite_element_problem/base_material.py +++ b/src/fenicsxconcrete/finite_element_problem/base_material.py @@ -92,16 +92,44 @@ def __init__( self.experiment = experiment self.mesh = self.experiment.mesh - # setting up default material parameters - default_fem_parameters = Parameters() - default_fem_parameters["g"] = 9.81 * ureg("m/s^2") - default_fem_parameters["dt"] = 1.0 * ureg("s") - - # adding experimental parameters to dictionary to combine to one - default_fem_parameters.update(self.experiment.parameters) + # initialize parameter attributes + setup_parameters = Parameters() + # setting up default setup parameters defined in each child + _, default_p = self.default_parameters() + setup_parameters.update(default_p) + # update with experiment parameters + setup_parameters.update(self.experiment.parameters) # update with input parameters - default_fem_parameters.update(parameters) - self.parameters = default_fem_parameters + setup_parameters.update(parameters) + + # get logger info which input parameters are set to default values + # plus check dimensionality of input parameters + keys_set_default = [] + for key in dict(default_p): + if key not in parameters: + keys_set_default.append(key) + else: + # check if units are compatible + dim_given = parameters[key].dimensionality + dim_default = default_p[key].dimensionality + print(f"check {key} {dim_given} {dim_default}") + if dim_given != dim_default: + raise ValueError( + f"given units for {key} are not compatible with default units: {dim_given} != {dim_default}" + ) + print(f"for the following parameters, the default values are used: {keys_set_default}") + self.logger.info(f"for the following parameters, the default values are used: {keys_set_default}") + + # # setting up default material parameters + # default_fem_parameters = Parameters() + # default_fem_parameters["g"] = 9.81 * ureg("m/s^2") + # default_fem_parameters["dt"] = 1.0 * ureg("s") + # + # # adding experimental parameters to dictionary to combine to one + # default_fem_parameters.update(self.experiment.parameters) + # # update with input parameters + # default_fem_parameters.update(parameters) + self.parameters = setup_parameters # remove units for use in fem model self.p = self.parameters.to_magnitude() self.experiment.p = self.p # update experimental parameter list for use in e.g. boundary definition @@ -135,6 +163,8 @@ def default_parameters() -> tuple[Experiment, dict[str, pint.Quantity]]: """returns a dictionary with required parameters and a set of working values as example""" # this must de defined in each setup class + pass + @abstractmethod def setup(self) -> None: # initialization of this specific problem diff --git a/src/fenicsxconcrete/finite_element_problem/concrete_am.py b/src/fenicsxconcrete/finite_element_problem/concrete_am.py index a168871..6c5a437 100644 --- a/src/fenicsxconcrete/finite_element_problem/concrete_am.py +++ b/src/fenicsxconcrete/finite_element_problem/concrete_am.py @@ -48,19 +48,19 @@ def __init__( """ - # adding default material parameter, will be overridden by outside input - default_p = Parameters() - default_p["stress_state"] = "plane_strain" * ureg("") # default stress state for 2D optional "plane_stress" - - # updating parameters, overriding defaults - default_p.update(parameters) + # # adding default material parameter, will be overridden by outside input + # default_p = Parameters() + # # default stress state for 2D optional "plane_stress" + # + # # updating parameters, overriding defaults + # default_p.update(parameters) if nonlinear_problem: self.nonlinear_problem = nonlinear_problem else: self.nonlinear_problem = ConcreteThixElasticModel # default material - super().__init__(experiment, default_p, pv_name, pv_path) + super().__init__(experiment, parameters, pv_name, pv_path) @staticmethod def parameter_description() -> dict[str, str]: @@ -73,12 +73,13 @@ def parameter_description() -> dict[str, str]: description = { "general parameters": { "rho": "density of fresh concrete", + "g": "gravity", "nu": "Poissons Ratio", "degree": "Polynomial degree for the FEM model", "q_degree": "Polynomial degree for which the quadrature rule integrates correctly", "load_time": "time in which the load is applied", "stress_state": "for 2D plain stress or plane strain", - "dt": "time step", # default set in material base class + "dt": "time step", }, "ThixElasticModel": { "E_0": "Youngs Modulus at age=0", @@ -112,10 +113,13 @@ def default_parameters( joined_parameters = { # Material parameter for concrete model with structural build-up "rho": 2070 * ureg("kg/m^3"), # density of fresh concrete + "g": 9.81 * ureg("m/s^2"), # gravity "nu": 0.3 * ureg(""), # Poissons Ratio # other model parameters - # "degree": 2 * ureg(""), # polynomial degree --> default defined in base_experiment!! + "degree": 2 * ureg(""), # polynomial degree "q_degree": 2 * ureg(""), # quadrature rule + "stress_state": "plane_strain" * ureg(""), # for 2D stress state + "dt": 1.0 * ureg("s"), # time step "load_time": 60 * ureg("s"), # body force load applied in s } if not non_linear_problem or non_linear_problem == ConcreteThixElasticModel: diff --git a/src/fenicsxconcrete/finite_element_problem/concrete_thermo_mechanical.py b/src/fenicsxconcrete/finite_element_problem/concrete_thermo_mechanical.py index 2f53779..eb5a238 100644 --- a/src/fenicsxconcrete/finite_element_problem/concrete_thermo_mechanical.py +++ b/src/fenicsxconcrete/finite_element_problem/concrete_thermo_mechanical.py @@ -32,20 +32,21 @@ def __init__( pv_path: str | None = None, ) -> None: - # adding default material parameter, will be overridden by outside input - default_p = Parameters() - # default_p['dummy'] = 'example' * ureg('') # example default parameter for this class + # # adding default material parameter, will be overridden by outside input + # default_p = Parameters() + # # default_p['dummy'] = 'example' * ureg('') # example default parameter for this class + # + # # updating parameters, overriding defaults + # default_p.update(parameters) - # updating parameters, overriding defaults - default_p.update(parameters) - - super().__init__(experiment, default_p, pv_name, pv_path) + super().__init__(experiment, parameters, pv_name, pv_path) @staticmethod def parameter_description() -> dict[str, str]: description = { "igc": "Ideal gas constant", "rho": "Density of concrete", + "g": "Gravitational acceleration", "themal_cond": "Thermal conductivity", "vol_heat_cap": "TODO", "Q_pot": "potential heat per weight of binder", @@ -68,6 +69,7 @@ def parameter_description() -> dict[str, str]: "ft_inf": "reference value for the tensile strength, default infinity, otherwise at alpha_tx", "a_ft": "exponential parameter to change the shape of the function ft(DOH)", "evolution_ft": "flag to turn off the evolution of the tensile strength", + "dt": "time step", } return description @@ -85,7 +87,8 @@ def default_parameters() -> tuple[Experiment, dict[str, pint.Quantity]]: default_parameters = { "igc": 8.3145 * ureg("J/K/mol"), "rho": 2350.0 * ureg("kg/m^3"), - "thermal_cond": 2.0 * ureg("W/(m*K)"), + "g": 9.81 * ureg("m/s^2"), + "thermal_cond": 2.0 * ureg("W/(m^3*K)"), "vol_heat_cap": 2.4e6 * ureg("J/(m^3 * K)"), # "Q_pot": 500e3 * ureg("J/kg"), only needed for postprocessing "Q_inf": 144000000 * ureg("J/m^3"), @@ -95,6 +98,7 @@ def default_parameters() -> tuple[Experiment, dict[str, pint.Quantity]]: "alpha_max": 0.875 * ureg(""), "alpha_tx": 0.68 * ureg(""), "T_ref": ureg.Quantity(25.0, ureg.degC), + "degree": 2 * ureg(""), "q_degree": 2 * ureg(""), "E_28": 15 * ureg("MPa"), "nu": 0.2 * ureg(""), @@ -106,6 +110,7 @@ def default_parameters() -> tuple[Experiment, dict[str, pint.Quantity]]: "ft_inf": 467000 * ureg(""), "a_ft": 1.0 * ureg(""), "evolution_ft": "True" * ureg(""), + "dt": 1.0 * ureg("s"), } default_parameters["E_act"] = 5653.0 * default_parameters["igc"] * ureg("J/mol") return experiment, default_parameters diff --git a/src/fenicsxconcrete/finite_element_problem/linear_elasticity.py b/src/fenicsxconcrete/finite_element_problem/linear_elasticity.py index bd32428..e2c3fef 100644 --- a/src/fenicsxconcrete/finite_element_problem/linear_elasticity.py +++ b/src/fenicsxconcrete/finite_element_problem/linear_elasticity.py @@ -21,14 +21,14 @@ def __init__( ) -> None: """defines default parameters, for the rest, see base class""" - # adding default material parameter, will be overridden by outside input - default_p = Parameters() - default_p["stress_state"] = "plane_strain" * ureg("") # default stress state in 2D, optional "plane_stress" + # # adding default material parameter, will be overridden by outside input + # default_p = Parameters() + # default_p["stress_state"] = "plane_strain" * ureg("") # default stress state in 2D, optional "plane_stress" + # + # # updating parameters, overriding defaults + # default_p.update(parameters) - # updating parameters, overriding defaults - default_p.update(parameters) - - super().__init__(experiment, default_p, pv_name, pv_path) + super().__init__(experiment, parameters, pv_name, pv_path) def setup(self) -> None: # compute different set of elastic moduli @@ -84,6 +84,27 @@ def setup(self) -> None: petsc_options={"ksp_type": "preonly", "pc_type": "lu"}, ) + @staticmethod + def parameter_description() -> dict[str, str]: + """static method returning a description dictionary for required parameters + + Returns: + description dictionary + + """ + description = { + "g": "gravity", + "dt": "time step", + "rho": "density of fresh concrete", + "E": "Young's Modulus", + "nu": "Poissons Ratio", + "stress_state": "for 2D plain stress or plane strain", + "degree": "Polynomial degree for the FEM model", + "dt": "time step", + } + + return description + @staticmethod def default_parameters() -> tuple[Experiment, dict[str, pint.Quantity]]: """returns a dictionary with required parameters and a set of working values as example""" @@ -91,10 +112,17 @@ def default_parameters() -> tuple[Experiment, dict[str, pint.Quantity]]: experiment = CantileverBeam(CantileverBeam.default_parameters()) model_parameters = {} + model_parameters["g"] = 9.81 * ureg("m/s^2") + model_parameters["dt"] = 1.0 * ureg("s") + model_parameters["rho"] = 7750 * ureg("kg/m^3") model_parameters["E"] = 210e9 * ureg("N/m^2") model_parameters["nu"] = 0.28 * ureg("") + model_parameters["stress_state"] = "plane_strain" * ureg("") + model_parameters["degree"] = 2 * ureg("") # polynomial degree + model_parameters["dt"] = 1.0 * ureg("s") + return experiment, model_parameters # Stress computation for linear elastic problem diff --git a/tests/finite_element_problem/test_linear_simple_cube.py b/tests/finite_element_problem/test_linear_simple_cube.py index 3b621f3..59d805c 100644 --- a/tests/finite_element_problem/test_linear_simple_cube.py +++ b/tests/finite_element_problem/test_linear_simple_cube.py @@ -153,3 +153,8 @@ def test_multiaxial_strain(dim: int, degree: int) -> None: assert result_corner == pytest.approx(target) assert result_center == pytest.approx(target / 2) + + +# main +if __name__ == "__main__": + test_disp(2) diff --git a/tests/finite_element_problem/test_thermo_mechanical_cube.py b/tests/finite_element_problem/test_thermo_mechanical_cube.py index 56e8da2..28944a2 100644 --- a/tests/finite_element_problem/test_thermo_mechanical_cube.py +++ b/tests/finite_element_problem/test_thermo_mechanical_cube.py @@ -159,7 +159,7 @@ def test_hydration_with_body_forces(dim: int): parameters["density_binder"] = 1440 * ureg("kg/m^3") # in kg/m^3 density of the binder parameters["thermal_cond"] = 2.0 * ureg( - "W/(m^3*K)" + "W/(m*K)" ) # effective thermal conductivity, approx in Wm^-3K^-1, concrete! # self.specific_heat_capacity = 9000 # effective specific heat capacity in J kg⁻1 K⁻1 parameters["vol_heat_cap"] = 2.4e6 * ureg("J/(m^3 * K)") # volumetric heat cap J/(m3 K)