diff --git a/tests/composite_viscous_outputs.cc b/tests/composite_viscous_outputs.cc index 542d2f29519..70a2f7c0c11 100644 --- a/tests/composite_viscous_outputs.cc +++ b/tests/composite_viscous_outputs.cc @@ -22,189 +22,193 @@ #include #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - // This function tests whether the composite creep rheology is producing - // the correct composite viscosity and partial strain rates corresponding to - // the different creep mechanisms incorporated into the rheology. - // It is assumed that each individual creep mechanism has already been tested. - - using namespace aspect::MaterialModel; - - // First, we set up a few objects which are used by the rheology model. - aspect::ParameterHandler prm; - const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); - auto n_phases = std::make_unique>(1); // 1 phase per composition - const unsigned int composition = 0; - const std::vector volume_fractions = {1.}; - const std::vector phase_function_values = std::vector(); - const std::vector n_phase_transitions_per_composition = std::vector(1); - - // Next, we initialise instances of the composite rheology and - // individual creep mechanisms. - std::unique_ptr> composite_creep; - composite_creep = std::make_unique>(); - composite_creep->initialize_simulator (simulator_access.get_simulator()); - composite_creep->declare_parameters(prm); - prm.set("Viscosity averaging scheme", "isostrain"); - prm.set("Include diffusion creep in composite rheology", "true"); - prm.set("Include dislocation creep in composite rheology", "true"); - prm.set("Include Peierls creep in composite rheology", "true"); - prm.set("Include Drucker Prager plasticity in composite rheology", "true"); - prm.set("Peierls creep flow law", "viscosity approximation"); - prm.set("Maximum yield stress", "5e8"); - composite_creep->parse_parameters(prm); - - std::unique_ptr> diffusion_creep; - diffusion_creep = std::make_unique>(); - diffusion_creep->initialize_simulator (simulator_access.get_simulator()); - diffusion_creep->declare_parameters(prm); - diffusion_creep->parse_parameters(prm); - - std::unique_ptr> dislocation_creep; - dislocation_creep = std::make_unique>(); - dislocation_creep->initialize_simulator (simulator_access.get_simulator()); - dislocation_creep->declare_parameters(prm); - dislocation_creep->parse_parameters(prm); - - std::unique_ptr> peierls_creep; - peierls_creep = std::make_unique>(); - peierls_creep->initialize_simulator (simulator_access.get_simulator()); - peierls_creep->declare_parameters(prm); - peierls_creep->parse_parameters(prm); - - std::unique_ptr> drucker_prager_power; - drucker_prager_power = std::make_unique>(); - drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); - drucker_prager_power->declare_parameters(prm); - prm.set("Maximum yield stress", "5e8"); - drucker_prager_power->parse_parameters(prm); - Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); - - // The creep components are arranged in series with each other. - // This package of components is then arranged in parallel with - // a strain rate limiter with a constant viscosity lim_visc. - // The whole system is then arranged in series with a viscosity limiter with - // viscosity max_visc. - // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) - double min_visc = prm.get_double("Minimum viscosity"); - double max_visc = prm.get_double("Maximum viscosity"); - double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); - - // Assign values to the variables which will be passed to compute_viscosity - // The test involves pure shear calculations at 1 GPa and variable temperature - double temperature; - const double pressure = 1.e9; - const double grain_size = 1.e-3; - SymmetricTensor<2,dim> strain_rate; - strain_rate[0][0] = -1e-11; - strain_rate[0][1] = 0.; - strain_rate[1][1] = 1e-11; - strain_rate[2][0] = 0.; - strain_rate[2][1] = 0.; - strain_rate[2][2] = 0.; - - std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; - - // Loop through strain rates, tracking whether there is a discrepancy in - // the decomposed strain rates. - bool error = false; - double viscosity; - double total_strain_rate; - double creep_strain_rate; - double creep_stress; - double diff_stress; - double disl_stress; - double prls_stress; - double drpr_stress; - std::vector partial_strain_rates(5, 0.); - - for (unsigned int i=0; i <= 10; i++) - { - temperature = 1000. + i*100.; - - // Compute the viscosity - viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); - total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); - - // The creep strain rate is calculated by subtracting the strain rate - // of the max viscosity dashpot from the total strain rate - // The creep stress is then calculated by subtracting the stress running - // through the strain rate limiter from the total stress - creep_strain_rate = total_strain_rate - partial_strain_rates[4]; - creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); - - // Print the output - std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; - for (unsigned int i=0; i < partial_strain_rates.size(); ++i) - { - std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; - } - std::cout << std::endl; - - // The following lines test that each individual creep mechanism - // experiences the same creep stress - - // Each creep mechanism should experience the same stress - diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); - disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); - prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); - if (partial_strain_rates[3] > 0.) - { - drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, - p.angle_internal_friction, - pressure, - partial_strain_rates[3], - p.max_yield_stress); - } - else - { - drpr_stress = creep_stress; - } - - if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) - { - error = true; - std::cout << " creep stress: " << creep_stress; - std::cout << " diffusion stress: " << diff_stress; - std::cout << " dislocation stress: " << disl_stress; - std::cout << " peierls stress: " << prls_stress; - std::cout << " drucker prager stress: " << drpr_stress << std::endl; - } - } - - if (error) - { - std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) + { + // This function tests whether the composite creep rheology is producing + // the correct composite viscosity and partial strain rates corresponding to + // the different creep mechanisms incorporated into the rheology. + // It is assumed that each individual creep mechanism has already been tested. + + using namespace aspect::MaterialModel; + + // First, we set up a few objects which are used by the rheology model. + aspect::ParameterHandler prm; + const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); + auto n_phases = std::make_unique>(1); // 1 phase per composition + const unsigned int composition = 0; + const std::vector volume_fractions = {1.}; + const std::vector phase_function_values = std::vector(); + const std::vector n_phase_transitions_per_composition = std::vector(1); + + // Next, we initialise instances of the composite rheology and + // individual creep mechanisms. + std::unique_ptr> composite_creep; + composite_creep = std::make_unique>(); + composite_creep->initialize_simulator (simulator_access.get_simulator()); + composite_creep->declare_parameters(prm); + prm.set("Viscosity averaging scheme", "isostrain"); + prm.set("Include diffusion creep in composite rheology", "true"); + prm.set("Include dislocation creep in composite rheology", "true"); + prm.set("Include Peierls creep in composite rheology", "true"); + prm.set("Include Drucker Prager plasticity in composite rheology", "true"); + prm.set("Peierls creep flow law", "viscosity approximation"); + prm.set("Maximum yield stress", "5e8"); + composite_creep->parse_parameters(prm); + + std::unique_ptr> diffusion_creep; + diffusion_creep = std::make_unique>(); + diffusion_creep->initialize_simulator (simulator_access.get_simulator()); + diffusion_creep->declare_parameters(prm); + diffusion_creep->parse_parameters(prm); + + std::unique_ptr> dislocation_creep; + dislocation_creep = std::make_unique>(); + dislocation_creep->initialize_simulator (simulator_access.get_simulator()); + dislocation_creep->declare_parameters(prm); + dislocation_creep->parse_parameters(prm); + + std::unique_ptr> peierls_creep; + peierls_creep = std::make_unique>(); + peierls_creep->initialize_simulator (simulator_access.get_simulator()); + peierls_creep->declare_parameters(prm); + peierls_creep->parse_parameters(prm); + + std::unique_ptr> drucker_prager_power; + drucker_prager_power = std::make_unique>(); + drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); + drucker_prager_power->declare_parameters(prm); + prm.set("Maximum yield stress", "5e8"); + drucker_prager_power->parse_parameters(prm); + Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + + // The creep components are arranged in series with each other. + // This package of components is then arranged in parallel with + // a strain rate limiter with a constant viscosity lim_visc. + // The whole system is then arranged in series with a viscosity limiter with + // viscosity max_visc. + // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) + double min_visc = prm.get_double("Minimum viscosity"); + double max_visc = prm.get_double("Maximum viscosity"); + double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); + + // Assign values to the variables which will be passed to compute_viscosity + // The test involves pure shear calculations at 1 GPa and variable temperature + double temperature; + const double pressure = 1.e9; + const double grain_size = 1.e-3; + SymmetricTensor<2,dim> strain_rate; + strain_rate[0][0] = -1e-11; + strain_rate[0][1] = 0.; + strain_rate[1][1] = 1e-11; + strain_rate[2][0] = 0.; + strain_rate[2][1] = 0.; + strain_rate[2][2] = 0.; + + std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; + + // Loop through strain rates, tracking whether there is a discrepancy in + // the decomposed strain rates. + bool error = false; + double viscosity; + double total_strain_rate; + double creep_strain_rate; + double creep_stress; + double diff_stress; + double disl_stress; + double prls_stress; + double drpr_stress; + std::vector partial_strain_rates(5, 0.); + + for (unsigned int i=0; i <= 10; i++) + { + temperature = 1000. + i*100.; + + // Compute the viscosity + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); + total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); + + // The creep strain rate is calculated by subtracting the strain rate + // of the max viscosity dashpot from the total strain rate + // The creep stress is then calculated by subtracting the stress running + // through the strain rate limiter from the total stress + creep_strain_rate = total_strain_rate - partial_strain_rates[4]; + creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); + + // Print the output + std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; + for (unsigned int i=0; i < partial_strain_rates.size(); ++i) + { + std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; + } + std::cout << std::endl; + + // The following lines test that each individual creep mechanism + // experiences the same creep stress + + // Each creep mechanism should experience the same stress + diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); + disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); + prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); + if (partial_strain_rates[3] > 0.) + { + drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, + p.angle_internal_friction, + pressure, + partial_strain_rates[3], + p.max_yield_stress); + } + else + { + drpr_stress = creep_stress; + } + + if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) + { + error = true; + std::cout << " creep stress: " << creep_stress; + std::cout << " diffusion stress: " << diff_stress; + std::cout << " dislocation stress: " << disl_stress; + std::cout << " peierls stress: " << prls_stress; + std::cout << " drucker prager stress: " << drpr_stress << std::endl; + } + } + + if (error) + { + std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + + } + + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } + + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } + + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &) -{ - AssertThrow(false,dealii::ExcInternalError()); -} - -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); -} - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/composite_viscous_outputs_isostress.cc b/tests/composite_viscous_outputs_isostress.cc index 110e3e861d9..dd6e30b0481 100644 --- a/tests/composite_viscous_outputs_isostress.cc +++ b/tests/composite_viscous_outputs_isostress.cc @@ -23,215 +23,219 @@ #include #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - // This function tests whether the composite creep rheology is producing - // the correct composite viscosity and partial strain rates corresponding to - // the different creep mechanisms incorporated into the rheology. - // It is assumed that each individual creep mechanism has already been tested. - - using namespace aspect::MaterialModel; - - // First, we set up a few objects which are used by the rheology model. - aspect::ParameterHandler prm; - - const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); - auto n_phases = std::make_unique>(1); // 1 phase per composition - const unsigned int composition = 0; - const std::vector volume_fractions = {0.6, 0.4}; - const std::vector phase_function_values = std::vector(); - const std::vector n_phase_transitions_per_composition = std::vector(1); - - // Next, we initialise instances of the composite rheology and - // individual creep mechanisms. - std::unique_ptr> composite_creep; - composite_creep = std::make_unique>(); - composite_creep->initialize_simulator (simulator_access.get_simulator()); - composite_creep->declare_parameters(prm); - prm.set("Viscosity averaging scheme", "isostress"); - prm.set("Include diffusion creep in composite rheology", "true"); - prm.set("Include dislocation creep in composite rheology", "true"); - prm.set("Include Peierls creep in composite rheology", "true"); - prm.set("Include Drucker Prager plasticity in composite rheology", "true"); - prm.set("Peierls creep flow law", "viscosity approximation"); - prm.set("Maximum yield stress", "5e8"); - composite_creep->parse_parameters(prm); - - std::unique_ptr> diffusion_creep; - diffusion_creep = std::make_unique>(); - diffusion_creep->initialize_simulator (simulator_access.get_simulator()); - diffusion_creep->declare_parameters(prm); - diffusion_creep->parse_parameters(prm); - - std::unique_ptr> dislocation_creep; - dislocation_creep = std::make_unique>(); - dislocation_creep->initialize_simulator (simulator_access.get_simulator()); - dislocation_creep->declare_parameters(prm); - dislocation_creep->parse_parameters(prm); - - std::unique_ptr> peierls_creep; - peierls_creep = std::make_unique>(); - peierls_creep->initialize_simulator (simulator_access.get_simulator()); - peierls_creep->declare_parameters(prm); - peierls_creep->parse_parameters(prm); - - std::unique_ptr> drucker_prager_power; - drucker_prager_power = std::make_unique>(); - drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); - drucker_prager_power->declare_parameters(prm); - prm.set("Maximum yield stress", "5e8"); - drucker_prager_power->parse_parameters(prm); - Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); - - // The creep components are arranged in series with each other. - // This package of components is then arranged in parallel with - // a strain rate limiter with a constant viscosity lim_visc. - // The whole system is then arranged in series with a viscosity limiter with - // viscosity max_visc. - // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) - double min_visc = prm.get_double("Minimum viscosity"); - double max_visc = prm.get_double("Maximum viscosity"); - double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); - - // Assign values to the variables which will be passed to compute_viscosity - // The test involves pure shear calculations at 1 GPa and variable temperature - double temperature; - const double pressure = 1.e9; - const double grain_size = 1.e-3; - SymmetricTensor<2,dim> strain_rate; - strain_rate[0][0] = -1e-11; - strain_rate[0][1] = 0.; - strain_rate[1][1] = 1e-11; - strain_rate[2][0] = 0.; - strain_rate[2][1] = 0.; - strain_rate[2][2] = 0.; - - std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; - - // Loop through strain rates, tracking whether there is a discrepancy in - // the decomposed strain rates. - bool error = false; - double viscosity; - double total_strain_rate; - double creep_strain_rate; - double creep_stress; - double diff_stress; - double disl_stress; - double prls_stress; - double drpr_stress; - std::vector partial_strain_rates(5, 0.); - - for (unsigned int i=0; i <= 10; i++) - { - temperature = 1000. + i*100.; - - // Compute the viscosity - viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); - total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); - - // The creep strain rate is calculated by subtracting the strain rate - // of the max viscosity dashpot from the total strain rate - // The creep stress is then calculated by subtracting the stress running - // through the strain rate limiter from the total stress - creep_strain_rate = total_strain_rate - partial_strain_rates[4]; - creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); - - // Print the output - std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; - for (unsigned int i=0; i < partial_strain_rates.size(); ++i) - { - std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; - } - std::cout << std::endl; - - // The following lines test that each individual creep mechanism - // experiences the same creep stress - - // Each creep mechanism should experience the same stress - diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); - disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); - prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); - if (partial_strain_rates[3] > 0.) - { - drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, - p.angle_internal_friction, - pressure, - partial_strain_rates[3], - p.max_yield_stress); - } - else - { - drpr_stress = creep_stress; - } - - if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) - { - error = true; - std::cout << " creep stress: " << creep_stress; - std::cout << " diffusion stress: " << diff_stress; - std::cout << " dislocation stress: " << disl_stress; - std::cout << " peierls stress: " << prls_stress; - std::cout << " drucker prager stress: " << drpr_stress << std::endl; - } - } - - if (error) - { - std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) + { + // This function tests whether the composite creep rheology is producing + // the correct composite viscosity and partial strain rates corresponding to + // the different creep mechanisms incorporated into the rheology. + // It is assumed that each individual creep mechanism has already been tested. + + using namespace aspect::MaterialModel; + + // First, we set up a few objects which are used by the rheology model. + aspect::ParameterHandler prm; + + const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); + auto n_phases = std::make_unique>(1); // 1 phase per composition + const unsigned int composition = 0; + const std::vector volume_fractions = {0.6, 0.4}; + const std::vector phase_function_values = std::vector(); + const std::vector n_phase_transitions_per_composition = std::vector(1); + + // Next, we initialise instances of the composite rheology and + // individual creep mechanisms. + std::unique_ptr> composite_creep; + composite_creep = std::make_unique>(); + composite_creep->initialize_simulator (simulator_access.get_simulator()); + composite_creep->declare_parameters(prm); + prm.set("Viscosity averaging scheme", "isostress"); + prm.set("Include diffusion creep in composite rheology", "true"); + prm.set("Include dislocation creep in composite rheology", "true"); + prm.set("Include Peierls creep in composite rheology", "true"); + prm.set("Include Drucker Prager plasticity in composite rheology", "true"); + prm.set("Peierls creep flow law", "viscosity approximation"); + prm.set("Maximum yield stress", "5e8"); + composite_creep->parse_parameters(prm); + + std::unique_ptr> diffusion_creep; + diffusion_creep = std::make_unique>(); + diffusion_creep->initialize_simulator (simulator_access.get_simulator()); + diffusion_creep->declare_parameters(prm); + diffusion_creep->parse_parameters(prm); + + std::unique_ptr> dislocation_creep; + dislocation_creep = std::make_unique>(); + dislocation_creep->initialize_simulator (simulator_access.get_simulator()); + dislocation_creep->declare_parameters(prm); + dislocation_creep->parse_parameters(prm); + + std::unique_ptr> peierls_creep; + peierls_creep = std::make_unique>(); + peierls_creep->initialize_simulator (simulator_access.get_simulator()); + peierls_creep->declare_parameters(prm); + peierls_creep->parse_parameters(prm); + + std::unique_ptr> drucker_prager_power; + drucker_prager_power = std::make_unique>(); + drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); + drucker_prager_power->declare_parameters(prm); + prm.set("Maximum yield stress", "5e8"); + drucker_prager_power->parse_parameters(prm); + Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + + // The creep components are arranged in series with each other. + // This package of components is then arranged in parallel with + // a strain rate limiter with a constant viscosity lim_visc. + // The whole system is then arranged in series with a viscosity limiter with + // viscosity max_visc. + // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) + double min_visc = prm.get_double("Minimum viscosity"); + double max_visc = prm.get_double("Maximum viscosity"); + double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); + + // Assign values to the variables which will be passed to compute_viscosity + // The test involves pure shear calculations at 1 GPa and variable temperature + double temperature; + const double pressure = 1.e9; + const double grain_size = 1.e-3; + SymmetricTensor<2,dim> strain_rate; + strain_rate[0][0] = -1e-11; + strain_rate[0][1] = 0.; + strain_rate[1][1] = 1e-11; + strain_rate[2][0] = 0.; + strain_rate[2][1] = 0.; + strain_rate[2][2] = 0.; + + std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; + + // Loop through strain rates, tracking whether there is a discrepancy in + // the decomposed strain rates. + bool error = false; + double viscosity; + double total_strain_rate; + double creep_strain_rate; + double creep_stress; + double diff_stress; + double disl_stress; + double prls_stress; + double drpr_stress; + std::vector partial_strain_rates(5, 0.); + + for (unsigned int i=0; i <= 10; i++) + { + temperature = 1000. + i*100.; + + // Compute the viscosity + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); + total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); + + // The creep strain rate is calculated by subtracting the strain rate + // of the max viscosity dashpot from the total strain rate + // The creep stress is then calculated by subtracting the stress running + // through the strain rate limiter from the total stress + creep_strain_rate = total_strain_rate - partial_strain_rates[4]; + creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); + + // Print the output + std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; + for (unsigned int i=0; i < partial_strain_rates.size(); ++i) + { + std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; + } + std::cout << std::endl; + + // The following lines test that each individual creep mechanism + // experiences the same creep stress + + // Each creep mechanism should experience the same stress + diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); + disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); + prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); + if (partial_strain_rates[3] > 0.) + { + drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, + p.angle_internal_friction, + pressure, + partial_strain_rates[3], + p.max_yield_stress); + } + else + { + drpr_stress = creep_stress; + } + + if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) + { + error = true; + std::cout << " creep stress: " << creep_stress; + std::cout << " diffusion stress: " << diff_stress; + std::cout << " dislocation stress: " << disl_stress; + std::cout << " peierls stress: " << prls_stress; + std::cout << " drucker prager stress: " << drpr_stress << std::endl; + } + } + + if (error) + { + std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } -} + } -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &) -{ - AssertThrow(false,dealii::ExcInternalError()); -} + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); -} + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } -using namespace aspect; + using namespace aspect; -void declare_parameters(const unsigned int dim, - ParameterHandler &prm) -{ - prm.enter_subsection("Compositional fields"); + void declare_parameters(const unsigned int dim, + ParameterHandler &prm) { - prm.declare_entry("Number of fields","1", Patterns::Integer()); + prm.enter_subsection("Compositional fields"); + { + prm.declare_entry("Number of fields","1", Patterns::Integer()); + } + prm.leave_subsection(); } - prm.leave_subsection(); -} -void parameter_connector () -{ - SimulatorSignals<2>::declare_additional_parameters.connect (&declare_parameters); - SimulatorSignals<3>::declare_additional_parameters.connect (&declare_parameters); -} + void parameter_connector () + { + SimulatorSignals<2>::declare_additional_parameters.connect (&declare_parameters); + SimulatorSignals<3>::declare_additional_parameters.connect (&declare_parameters); + } + + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) + ASPECT_REGISTER_SIGNALS_PARAMETER_CONNECTOR(parameter_connector) -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) -ASPECT_REGISTER_SIGNALS_PARAMETER_CONNECTOR(parameter_connector) +} diff --git a/tests/composite_viscous_outputs_limited.cc b/tests/composite_viscous_outputs_limited.cc index 08bbf4b895d..fdcbf3f551d 100644 --- a/tests/composite_viscous_outputs_limited.cc +++ b/tests/composite_viscous_outputs_limited.cc @@ -22,200 +22,204 @@ #include #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - // This function tests whether the composite creep rheology is producing - // the correct composite viscosity and partial strain rates corresponding to - // the different creep mechanisms incorporated into the rheology. - // It is assumed that each individual creep mechanism has already been tested. - - using namespace aspect::MaterialModel; - - // First, we set up a few objects which are used by the rheology model. - aspect::ParameterHandler prm; - const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); - auto n_phases = std::make_unique>(1); // 1 phase per composition - const unsigned int composition = 0; - const std::vector volume_fractions = {1.}; - const std::vector phase_function_values = std::vector(); - const std::vector n_phase_transitions_per_composition = std::vector(1); - - // Next, we initialise instances of the composite rheology and - // individual creep mechanisms. - std::unique_ptr> composite_creep; - composite_creep = std::make_unique>(); - composite_creep->initialize_simulator (simulator_access.get_simulator()); - composite_creep->declare_parameters(prm); - prm.set("Viscosity averaging scheme", "isostrain"); - prm.set("Include diffusion creep in composite rheology", "true"); - prm.set("Include dislocation creep in composite rheology", "true"); - prm.set("Include Peierls creep in composite rheology", "true"); - prm.set("Include Drucker Prager plasticity in composite rheology", "true"); - prm.set("Peierls creep flow law", "viscosity approximation"); - prm.set("Maximum yield stress", "5e8"); - prm.set("Minimum viscosity", "1e18"); - prm.set("Maximum viscosity", "4e18"); - composite_creep->parse_parameters(prm); - - std::unique_ptr> diffusion_creep; - diffusion_creep = std::make_unique>(); - diffusion_creep->initialize_simulator (simulator_access.get_simulator()); - diffusion_creep->declare_parameters(prm); - diffusion_creep->parse_parameters(prm); - - std::unique_ptr> dislocation_creep; - dislocation_creep = std::make_unique>(); - dislocation_creep->initialize_simulator (simulator_access.get_simulator()); - dislocation_creep->declare_parameters(prm); - dislocation_creep->parse_parameters(prm); - - std::unique_ptr> peierls_creep; - peierls_creep = std::make_unique>(); - peierls_creep->initialize_simulator (simulator_access.get_simulator()); - peierls_creep->declare_parameters(prm); - peierls_creep->parse_parameters(prm); - - std::unique_ptr> drucker_prager_power; - drucker_prager_power = std::make_unique>(); - drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); - drucker_prager_power->declare_parameters(prm); - prm.set("Maximum yield stress", "5e8"); - drucker_prager_power->parse_parameters(prm); - Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); - - // The creep components are arranged in series with each other. - // This package of components is then arranged in parallel with - // a strain rate limiter with a constant viscosity lim_visc. - // The whole system is then arranged in series with a viscosity limiter with - // viscosity maximum_viscosity. - // lim_visc is equal to (minimum_viscosity*maximum_viscosity)/(maximum_viscosity - minimum_viscosity) - double minimum_viscosity = prm.get_double("Minimum viscosity"); - double maximum_viscosity = prm.get_double("Maximum viscosity"); - - AssertThrow(minimum_viscosity > 0, - ExcMessage("Minimum viscosity needs to be larger than zero.")); - - AssertThrow(maximum_viscosity > 1.1 * minimum_viscosity, - ExcMessage("The maximum viscosity needs to be at least ten percent larger than the minimum viscosity. " - "If you require an isoviscous model consider a different rheology, or set the " - "parameters of the active flow laws to be independent of temperature, pressure, grain size, and stress.")); - - double lim_visc = (minimum_viscosity*maximum_viscosity)/(maximum_viscosity - minimum_viscosity); - - // Assign values to the variables which will be passed to compute_viscosity - // The test involves pure shear calculations at 1 GPa and variable temperature - double temperature; - const double pressure = 1.e9; - const double grain_size = 1.e-3; - SymmetricTensor<2,dim> strain_rate; - strain_rate[0][0] = -1e-11; - strain_rate[0][1] = 0.; - strain_rate[1][1] = 1e-11; - strain_rate[2][0] = 0.; - strain_rate[2][1] = 0.; - strain_rate[2][2] = 0.; - - std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; - - // Loop through strain rates, tracking whether there is a discrepancy in - // the decomposed strain rates. - bool error = false; - double viscosity; - double total_strain_rate; - double creep_strain_rate; - double creep_stress; - double diff_stress; - double disl_stress; - double prls_stress; - double drpr_stress; - std::vector partial_strain_rates(5, 0.); - - for (unsigned int i=0; i <= 10; i++) - { - temperature = 1000. + i*100.; - - // Compute the viscosity - viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); - total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); - - // The creep strain rate is calculated by subtracting the strain rate - // of the max viscosity dashpot from the total strain rate - // The creep stress is then calculated by subtracting the stress running - // through the strain rate limiter from the total stress - creep_strain_rate = total_strain_rate - partial_strain_rates[4]; - creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); - - // Print the output - std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; - for (unsigned int i=0; i < partial_strain_rates.size(); ++i) - { - std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; - } - std::cout << std::endl; - - // The following lines test that each individual creep mechanism - // experiences the same creep stress - - // Each creep mechanism should experience the same stress - diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); - disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); - prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); - if (partial_strain_rates[3] > 0.) - { - drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, - p.angle_internal_friction, - pressure, - partial_strain_rates[3], - p.max_yield_stress); - } - else - { - drpr_stress = creep_stress; - } - - if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) - { - error = true; - std::cout << " creep stress: " << creep_stress; - std::cout << " diffusion stress: " << diff_stress; - std::cout << " dislocation stress: " << disl_stress; - std::cout << " peierls stress: " << prls_stress; - std::cout << " drucker prager stress: " << drpr_stress << std::endl; - } - } - - if (error) - { - std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) + { + // This function tests whether the composite creep rheology is producing + // the correct composite viscosity and partial strain rates corresponding to + // the different creep mechanisms incorporated into the rheology. + // It is assumed that each individual creep mechanism has already been tested. + + using namespace aspect::MaterialModel; + + // First, we set up a few objects which are used by the rheology model. + aspect::ParameterHandler prm; + const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); + auto n_phases = std::make_unique>(1); // 1 phase per composition + const unsigned int composition = 0; + const std::vector volume_fractions = {1.}; + const std::vector phase_function_values = std::vector(); + const std::vector n_phase_transitions_per_composition = std::vector(1); + + // Next, we initialise instances of the composite rheology and + // individual creep mechanisms. + std::unique_ptr> composite_creep; + composite_creep = std::make_unique>(); + composite_creep->initialize_simulator (simulator_access.get_simulator()); + composite_creep->declare_parameters(prm); + prm.set("Viscosity averaging scheme", "isostrain"); + prm.set("Include diffusion creep in composite rheology", "true"); + prm.set("Include dislocation creep in composite rheology", "true"); + prm.set("Include Peierls creep in composite rheology", "true"); + prm.set("Include Drucker Prager plasticity in composite rheology", "true"); + prm.set("Peierls creep flow law", "viscosity approximation"); + prm.set("Maximum yield stress", "5e8"); + prm.set("Minimum viscosity", "1e18"); + prm.set("Maximum viscosity", "4e18"); + composite_creep->parse_parameters(prm); + + std::unique_ptr> diffusion_creep; + diffusion_creep = std::make_unique>(); + diffusion_creep->initialize_simulator (simulator_access.get_simulator()); + diffusion_creep->declare_parameters(prm); + diffusion_creep->parse_parameters(prm); + + std::unique_ptr> dislocation_creep; + dislocation_creep = std::make_unique>(); + dislocation_creep->initialize_simulator (simulator_access.get_simulator()); + dislocation_creep->declare_parameters(prm); + dislocation_creep->parse_parameters(prm); + + std::unique_ptr> peierls_creep; + peierls_creep = std::make_unique>(); + peierls_creep->initialize_simulator (simulator_access.get_simulator()); + peierls_creep->declare_parameters(prm); + peierls_creep->parse_parameters(prm); + + std::unique_ptr> drucker_prager_power; + drucker_prager_power = std::make_unique>(); + drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); + drucker_prager_power->declare_parameters(prm); + prm.set("Maximum yield stress", "5e8"); + drucker_prager_power->parse_parameters(prm); + Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + + // The creep components are arranged in series with each other. + // This package of components is then arranged in parallel with + // a strain rate limiter with a constant viscosity lim_visc. + // The whole system is then arranged in series with a viscosity limiter with + // viscosity maximum_viscosity. + // lim_visc is equal to (minimum_viscosity*maximum_viscosity)/(maximum_viscosity - minimum_viscosity) + double minimum_viscosity = prm.get_double("Minimum viscosity"); + double maximum_viscosity = prm.get_double("Maximum viscosity"); + + AssertThrow(minimum_viscosity > 0, + ExcMessage("Minimum viscosity needs to be larger than zero.")); + + AssertThrow(maximum_viscosity > 1.1 * minimum_viscosity, + ExcMessage("The maximum viscosity needs to be at least ten percent larger than the minimum viscosity. " + "If you require an isoviscous model consider a different rheology, or set the " + "parameters of the active flow laws to be independent of temperature, pressure, grain size, and stress.")); + + double lim_visc = (minimum_viscosity*maximum_viscosity)/(maximum_viscosity - minimum_viscosity); + + // Assign values to the variables which will be passed to compute_viscosity + // The test involves pure shear calculations at 1 GPa and variable temperature + double temperature; + const double pressure = 1.e9; + const double grain_size = 1.e-3; + SymmetricTensor<2,dim> strain_rate; + strain_rate[0][0] = -1e-11; + strain_rate[0][1] = 0.; + strain_rate[1][1] = 1e-11; + strain_rate[2][0] = 0.; + strain_rate[2][1] = 0.; + strain_rate[2][2] = 0.; + + std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; + + // Loop through strain rates, tracking whether there is a discrepancy in + // the decomposed strain rates. + bool error = false; + double viscosity; + double total_strain_rate; + double creep_strain_rate; + double creep_stress; + double diff_stress; + double disl_stress; + double prls_stress; + double drpr_stress; + std::vector partial_strain_rates(5, 0.); + + for (unsigned int i=0; i <= 10; i++) + { + temperature = 1000. + i*100.; + + // Compute the viscosity + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); + total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); + + // The creep strain rate is calculated by subtracting the strain rate + // of the max viscosity dashpot from the total strain rate + // The creep stress is then calculated by subtracting the stress running + // through the strain rate limiter from the total stress + creep_strain_rate = total_strain_rate - partial_strain_rates[4]; + creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); + + // Print the output + std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; + for (unsigned int i=0; i < partial_strain_rates.size(); ++i) + { + std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; + } + std::cout << std::endl; + + // The following lines test that each individual creep mechanism + // experiences the same creep stress + + // Each creep mechanism should experience the same stress + diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); + disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); + prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); + if (partial_strain_rates[3] > 0.) + { + drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, + p.angle_internal_friction, + pressure, + partial_strain_rates[3], + p.max_yield_stress); + } + else + { + drpr_stress = creep_stress; + } + + if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) + { + error = true; + std::cout << " creep stress: " << creep_stress; + std::cout << " diffusion stress: " << diff_stress; + std::cout << " dislocation stress: " << disl_stress; + std::cout << " peierls stress: " << prls_stress; + std::cout << " drucker prager stress: " << drpr_stress << std::endl; + } + } + + if (error) + { + std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + + } + + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } + + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } + + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &) -{ - AssertThrow(false,dealii::ExcInternalError()); -} - -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); -} - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/composite_viscous_outputs_no_peierls.cc b/tests/composite_viscous_outputs_no_peierls.cc index ef6ef4c702b..c7826412354 100644 --- a/tests/composite_viscous_outputs_no_peierls.cc +++ b/tests/composite_viscous_outputs_no_peierls.cc @@ -21,180 +21,183 @@ #include #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - // This function tests whether the composite creep rheology is producing - // the correct composite viscosity and partial strain rates corresponding to - // the different creep mechanisms incorporated into the rheology. - // It is assumed that each individual creep mechanism has already been tested. - - using namespace aspect::MaterialModel; - - // First, we set up a few objects which are used by the rheology model. - aspect::ParameterHandler prm; - - const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); - auto n_phases = std::make_unique>(1); // 1 phase per composition - const unsigned int composition = 0; - const std::vector volume_fractions = {1.}; - const std::vector phase_function_values = std::vector(); - const std::vector n_phase_transitions_per_composition = std::vector(1); - - // Next, we initialise instances of the composite rheology and - // individual creep mechanisms. - std::unique_ptr> composite_creep; - composite_creep = std::make_unique>(); - composite_creep->initialize_simulator (simulator_access.get_simulator()); - composite_creep->declare_parameters(prm); - prm.set("Viscosity averaging scheme", "isostrain"); - prm.set("Include diffusion creep in composite rheology", "true"); - prm.set("Include dislocation creep in composite rheology", "true"); - prm.set("Include Peierls creep in composite rheology", "false"); - prm.set("Include Drucker Prager plasticity in composite rheology", "true"); - prm.set("Peierls creep flow law", "viscosity approximation"); - prm.set("Maximum yield stress", "5e8"); - composite_creep->parse_parameters(prm); - - std::unique_ptr> diffusion_creep; - diffusion_creep = std::make_unique>(); - diffusion_creep->initialize_simulator (simulator_access.get_simulator()); - diffusion_creep->declare_parameters(prm); - diffusion_creep->parse_parameters(prm); - - std::unique_ptr> dislocation_creep; - dislocation_creep = std::make_unique>(); - dislocation_creep->initialize_simulator (simulator_access.get_simulator()); - dislocation_creep->declare_parameters(prm); - dislocation_creep->parse_parameters(prm); - - std::unique_ptr> drucker_prager_power; - drucker_prager_power = std::make_unique>(); - drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); - drucker_prager_power->declare_parameters(prm); - prm.set("Maximum yield stress", "5e8"); - drucker_prager_power->parse_parameters(prm); - Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); - - // The creep components are arranged in series with each other. - // This package of components is then arranged in parallel with - // a strain rate limiter with a constant viscosity lim_visc. - // The whole system is then arranged in series with a viscosity limiter with - // viscosity max_visc. - // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) - double min_visc = prm.get_double("Minimum viscosity"); - double max_visc = prm.get_double("Maximum viscosity"); - double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); - - // Assign values to the variables which will be passed to compute_viscosity - // The test involves pure shear calculations at 1 GPa and variable temperature - double temperature; - const double pressure = 1.e9; - const double grain_size = 1.e-3; - SymmetricTensor<2,dim> strain_rate; - strain_rate[0][0] = -1e-11; - strain_rate[0][1] = 0.; - strain_rate[1][1] = 1e-11; - strain_rate[2][0] = 0.; - strain_rate[2][1] = 0.; - strain_rate[2][2] = 0.; - - std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; - - // Loop through strain rates, tracking whether there is a discrepancy in - // the decomposed strain rates. - bool error = false; - double viscosity; - double total_strain_rate; - double creep_strain_rate; - double creep_stress; - double diff_stress; - double disl_stress; - double drpr_stress; - std::vector partial_strain_rates(5, 0.); - - for (unsigned int i=0; i <= 10; i++) - { - temperature = 1000. + i*100.; - - // Compute the viscosity - viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); - total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); - - // The creep strain rate is calculated by subtracting the strain rate - // of the max viscosity dashpot from the total strain rate - // The creep stress is then calculated by subtracting the stress running - // through the strain rate limiter from the total stress - creep_strain_rate = total_strain_rate - partial_strain_rates[4]; - creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); - - // Print the output - std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; - for (unsigned int i=0; i < partial_strain_rates.size(); ++i) - { - std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; - } - std::cout << std::endl; - - // The following lines test that each individual creep mechanism - // experiences the same creep stress - - // Each creep mechanism should experience the same stress - diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); - disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); - if (partial_strain_rates[3] > 0.) - { - drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, - p.angle_internal_friction, - pressure, - partial_strain_rates[3], - p.max_yield_stress); - } - else - { - drpr_stress = creep_stress; - } - - if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) - { - error = true; - std::cout << " creep stress: " << creep_stress; - std::cout << " diffusion stress: " << diff_stress; - std::cout << " dislocation stress: " << disl_stress; - std::cout << " drucker prager stress: " << drpr_stress << std::endl; - } - } - - if (error) - { - std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } - -} - -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &) -{ - AssertThrow(false,dealii::ExcInternalError()); + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) + { + // This function tests whether the composite creep rheology is producing + // the correct composite viscosity and partial strain rates corresponding to + // the different creep mechanisms incorporated into the rheology. + // It is assumed that each individual creep mechanism has already been tested. + + using namespace aspect::MaterialModel; + + // First, we set up a few objects which are used by the rheology model. + aspect::ParameterHandler prm; + + const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); + auto n_phases = std::make_unique>(1); // 1 phase per composition + const unsigned int composition = 0; + const std::vector volume_fractions = {1.}; + const std::vector phase_function_values = std::vector(); + const std::vector n_phase_transitions_per_composition = std::vector(1); + + // Next, we initialise instances of the composite rheology and + // individual creep mechanisms. + std::unique_ptr> composite_creep; + composite_creep = std::make_unique>(); + composite_creep->initialize_simulator (simulator_access.get_simulator()); + composite_creep->declare_parameters(prm); + prm.set("Viscosity averaging scheme", "isostrain"); + prm.set("Include diffusion creep in composite rheology", "true"); + prm.set("Include dislocation creep in composite rheology", "true"); + prm.set("Include Peierls creep in composite rheology", "false"); + prm.set("Include Drucker Prager plasticity in composite rheology", "true"); + prm.set("Peierls creep flow law", "viscosity approximation"); + prm.set("Maximum yield stress", "5e8"); + composite_creep->parse_parameters(prm); + + std::unique_ptr> diffusion_creep; + diffusion_creep = std::make_unique>(); + diffusion_creep->initialize_simulator (simulator_access.get_simulator()); + diffusion_creep->declare_parameters(prm); + diffusion_creep->parse_parameters(prm); + + std::unique_ptr> dislocation_creep; + dislocation_creep = std::make_unique>(); + dislocation_creep->initialize_simulator (simulator_access.get_simulator()); + dislocation_creep->declare_parameters(prm); + dislocation_creep->parse_parameters(prm); + + std::unique_ptr> drucker_prager_power; + drucker_prager_power = std::make_unique>(); + drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); + drucker_prager_power->declare_parameters(prm); + prm.set("Maximum yield stress", "5e8"); + drucker_prager_power->parse_parameters(prm); + Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + + // The creep components are arranged in series with each other. + // This package of components is then arranged in parallel with + // a strain rate limiter with a constant viscosity lim_visc. + // The whole system is then arranged in series with a viscosity limiter with + // viscosity max_visc. + // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) + double min_visc = prm.get_double("Minimum viscosity"); + double max_visc = prm.get_double("Maximum viscosity"); + double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); + + // Assign values to the variables which will be passed to compute_viscosity + // The test involves pure shear calculations at 1 GPa and variable temperature + double temperature; + const double pressure = 1.e9; + const double grain_size = 1.e-3; + SymmetricTensor<2,dim> strain_rate; + strain_rate[0][0] = -1e-11; + strain_rate[0][1] = 0.; + strain_rate[1][1] = 1e-11; + strain_rate[2][0] = 0.; + strain_rate[2][1] = 0.; + strain_rate[2][2] = 0.; + + std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; + + // Loop through strain rates, tracking whether there is a discrepancy in + // the decomposed strain rates. + bool error = false; + double viscosity; + double total_strain_rate; + double creep_strain_rate; + double creep_stress; + double diff_stress; + double disl_stress; + double drpr_stress; + std::vector partial_strain_rates(5, 0.); + + for (unsigned int i=0; i <= 10; i++) + { + temperature = 1000. + i*100.; + + // Compute the viscosity + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates); + total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); + + // The creep strain rate is calculated by subtracting the strain rate + // of the max viscosity dashpot from the total strain rate + // The creep stress is then calculated by subtracting the stress running + // through the strain rate limiter from the total stress + creep_strain_rate = total_strain_rate - partial_strain_rates[4]; + creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); + + // Print the output + std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; + for (unsigned int i=0; i < partial_strain_rates.size(); ++i) + { + std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; + } + std::cout << std::endl; + + // The following lines test that each individual creep mechanism + // experiences the same creep stress + + // Each creep mechanism should experience the same stress + diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); + disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); + if (partial_strain_rates[3] > 0.) + { + drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, + p.angle_internal_friction, + pressure, + partial_strain_rates[3], + p.max_yield_stress); + } + else + { + drpr_stress = creep_stress; + } + + if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) + { + error = true; + std::cout << " creep stress: " << creep_stress; + std::cout << " diffusion stress: " << diff_stress; + std::cout << " dislocation stress: " << disl_stress; + std::cout << " drucker prager stress: " << drpr_stress << std::endl; + } + } + + if (error) + { + std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + + } + + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } + + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } + + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); -} - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/composite_viscous_outputs_phases.cc b/tests/composite_viscous_outputs_phases.cc index c4fc9e024c9..932bb762ed7 100644 --- a/tests/composite_viscous_outputs_phases.cc +++ b/tests/composite_viscous_outputs_phases.cc @@ -23,238 +23,241 @@ #include #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - // This function tests whether the composite creep rheology is producing - // the correct composite viscosity and partial strain rates corresponding to - // the different creep mechanisms incorporated into the rheology. - // It is assumed that each individual creep mechanism has already been tested. - - using namespace aspect::MaterialModel; - - // First, we set up a few objects which are used by the rheology model. - aspect::ParameterHandler prm; - const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); - MaterialUtilities::PhaseFunction phase_function; - phase_function.initialize_simulator (simulator_access.get_simulator()); - phase_function.declare_parameters (prm); - prm.set("Define transition by depth instead of pressure", "false"); - prm.set("Phase transition pressures", "3e9"); - prm.set("Phase transition pressure widths", "1e9"); - prm.set("Phase transition temperatures", "273"); - prm.set("Phase transition Clapeyron slopes", "0"); - phase_function.parse_parameters (prm); - - std::vector n_phases_for_each_composition = phase_function.n_phases_for_each_composition(); - // Currently, phase_function.n_phases_for_each_composition() returns a list of length - // equal to the total number of compositions, whether or not they are chemical compositions. - // The equation_of_state (multicomponent incompressible) requires a list only for - // chemical compositions. - std::vector n_phases_for_each_chemical_composition = {n_phases_for_each_composition[0]}; - std::vector n_phase_transitions_for_each_chemical_composition = {n_phases_for_each_composition[0] - 1}; - unsigned int n_phases = n_phases_for_each_composition[0]; - for (auto i : simulator_access.introspection().chemical_composition_field_indices()) - { - n_phases_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1]); - n_phase_transitions_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1] - 1); - n_phases += n_phases_for_each_composition[i+1]; - } - - const unsigned int composition = 0; - const std::vector volume_fractions = {1.}; - std::vector phase_function_values = {0.}; - const std::vector n_phase_transitions_per_composition = n_phase_transitions_for_each_chemical_composition; - - // Next, we initialise instances of the composite rheology and - // individual creep mechanisms. - std::unique_ptr> composite_creep; - composite_creep = std::make_unique>(); - composite_creep->initialize_simulator (simulator_access.get_simulator()); - composite_creep->declare_parameters(prm); - MaterialUtilities::PhaseFunction::declare_parameters(prm); - prm.set("Viscosity averaging scheme", "isostrain"); - prm.set("Include diffusion creep in composite rheology", "true"); - prm.set("Include dislocation creep in composite rheology", "true"); - prm.set("Include Peierls creep in composite rheology", "true"); - prm.set("Include Drucker Prager plasticity in composite rheology", "true"); - prm.set("Peierls creep flow law", "viscosity approximation"); - prm.set("Cohesions", "background:1e9|5e8"); - prm.set("Maximum yield stress", "5e10"); - composite_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); - - std::unique_ptr> diffusion_creep; - diffusion_creep = std::make_unique>(); - diffusion_creep->initialize_simulator (simulator_access.get_simulator()); - diffusion_creep->declare_parameters(prm); - diffusion_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); - - std::unique_ptr> dislocation_creep; - dislocation_creep = std::make_unique>(); - dislocation_creep->initialize_simulator (simulator_access.get_simulator()); - dislocation_creep->declare_parameters(prm); - dislocation_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); - - std::unique_ptr> peierls_creep; - peierls_creep = std::make_unique>(); - peierls_creep->initialize_simulator (simulator_access.get_simulator()); - peierls_creep->declare_parameters(prm); - peierls_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); - - std::unique_ptr> drucker_prager_power; - drucker_prager_power = std::make_unique>(); - drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); - drucker_prager_power->declare_parameters(prm); - prm.set("Cohesions", "background:1e9|5e8"); - prm.set("Maximum yield stress", "5e10"); - drucker_prager_power->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); - - // The creep components are arranged in series with each other. - // This package of components is then arranged in parallel with - // a strain rate limiter with a constant viscosity lim_visc. - // The whole system is then arranged in series with a viscosity limiter with - // viscosity max_visc. - // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) - double min_visc = prm.get_double("Minimum viscosity"); - double max_visc = prm.get_double("Maximum viscosity"); - double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); - - // Assign values to the variables which will be passed to compute_viscosity - // The test involves pure shear calculations at variable pressure and temperature - double temperature; - double pressure; - const double grain_size = 1.e-3; - SymmetricTensor<2,dim> strain_rate; - strain_rate[0][0] = -1e-11; - strain_rate[0][1] = 0.; - strain_rate[1][1] = 1e-11; - strain_rate[2][0] = 0.; - strain_rate[2][1] = 0.; - strain_rate[2][2] = 0.; - - std::cout << "temperature (K) phase transition progress eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; - - // Loop through strain rates, tracking whether there is a discrepancy in - // the decomposed strain rates. - bool error = false; - double viscosity; - double total_strain_rate; - double creep_strain_rate; - double creep_stress; - double diff_stress; - double disl_stress; - double prls_stress; - double drpr_stress; - std::vector partial_strain_rates(5, 0.); - - for (unsigned int i=0; i <= 2; i++) - { - pressure = 1.e9 + i*2.e9; - std::cout << "pressure: " << pressure / 1.e9 << " GPa" << std::endl; - for (unsigned int j = 0; j <= 10; j++) - { - temperature = 1000. + j*100.; - - // Compute the phase function values - // The depth and gravity are set to zero because they are unused - // when phase functions are calculated by pressure. - // The phase index is set to invalid_unsigned_int, because it is only used internally - // in phase_average_equation_of_state_outputs to loop over all existing phases - MaterialUtilities::PhaseFunctionInputs phase_inputs(temperature, - pressure, - 0., 0., - numbers::invalid_unsigned_int); - - // Compute value of phase functions - for (unsigned int j=0; j < phase_function.n_phase_transitions(); ++j) - { - phase_inputs.phase_transition_index = j; - phase_function_values[j] = phase_function.compute_value(phase_inputs); - } - - Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); - - // Compute the viscosity - viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates, phase_function_values, n_phase_transitions_per_composition); - total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); - - // The creep strain rate is calculated by subtracting the strain rate - // of the max viscosity dashpot from the total strain rate - // The creep stress is then calculated by subtracting the stress running - // through the strain rate limiter from the total stress - creep_strain_rate = total_strain_rate - partial_strain_rates[4]; - creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); - - // Print the output - std::cout << temperature << ' ' << phase_function_values[0] << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; - for (unsigned int i=0; i < partial_strain_rates.size(); ++i) - { - std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; - } - std::cout << std::endl; - - // The following lines test that each individual creep mechanism - // experiences the same creep stress - - // Each creep mechanism should experience the same stress - diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); - disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); - prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); - if (partial_strain_rates[3] > 0.) - { - drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, - p.angle_internal_friction, - pressure, - partial_strain_rates[3], - p.max_yield_stress); - } - else - { - drpr_stress = creep_stress; - } - - if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) - || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) - { - error = true; - std::cout << " creep stress: " << creep_stress; - std::cout << " diffusion stress: " << diff_stress; - std::cout << " dislocation stress: " << disl_stress; - std::cout << " peierls stress: " << prls_stress; - std::cout << " drucker prager stress: " << drpr_stress << std::endl; - } - } - } - - if (error) - { - std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) + { + // This function tests whether the composite creep rheology is producing + // the correct composite viscosity and partial strain rates corresponding to + // the different creep mechanisms incorporated into the rheology. + // It is assumed that each individual creep mechanism has already been tested. + + using namespace aspect::MaterialModel; + + // First, we set up a few objects which are used by the rheology model. + aspect::ParameterHandler prm; + const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names(); + MaterialUtilities::PhaseFunction phase_function; + phase_function.initialize_simulator (simulator_access.get_simulator()); + phase_function.declare_parameters (prm); + prm.set("Define transition by depth instead of pressure", "false"); + prm.set("Phase transition pressures", "3e9"); + prm.set("Phase transition pressure widths", "1e9"); + prm.set("Phase transition temperatures", "273"); + prm.set("Phase transition Clapeyron slopes", "0"); + phase_function.parse_parameters (prm); + + std::vector n_phases_for_each_composition = phase_function.n_phases_for_each_composition(); + // Currently, phase_function.n_phases_for_each_composition() returns a list of length + // equal to the total number of compositions, whether or not they are chemical compositions. + // The equation_of_state (multicomponent incompressible) requires a list only for + // chemical compositions. + std::vector n_phases_for_each_chemical_composition = {n_phases_for_each_composition[0]}; + std::vector n_phase_transitions_for_each_chemical_composition = {n_phases_for_each_composition[0] - 1}; + unsigned int n_phases = n_phases_for_each_composition[0]; + for (auto i : simulator_access.introspection().chemical_composition_field_indices()) + { + n_phases_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1]); + n_phase_transitions_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1] - 1); + n_phases += n_phases_for_each_composition[i+1]; + } + + const unsigned int composition = 0; + const std::vector volume_fractions = {1.}; + std::vector phase_function_values = {0.}; + const std::vector n_phase_transitions_per_composition = n_phase_transitions_for_each_chemical_composition; + + // Next, we initialise instances of the composite rheology and + // individual creep mechanisms. + std::unique_ptr> composite_creep; + composite_creep = std::make_unique>(); + composite_creep->initialize_simulator (simulator_access.get_simulator()); + composite_creep->declare_parameters(prm); + MaterialUtilities::PhaseFunction::declare_parameters(prm); + prm.set("Viscosity averaging scheme", "isostrain"); + prm.set("Include diffusion creep in composite rheology", "true"); + prm.set("Include dislocation creep in composite rheology", "true"); + prm.set("Include Peierls creep in composite rheology", "true"); + prm.set("Include Drucker Prager plasticity in composite rheology", "true"); + prm.set("Peierls creep flow law", "viscosity approximation"); + prm.set("Cohesions", "background:1e9|5e8"); + prm.set("Maximum yield stress", "5e10"); + composite_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); + + std::unique_ptr> diffusion_creep; + diffusion_creep = std::make_unique>(); + diffusion_creep->initialize_simulator (simulator_access.get_simulator()); + diffusion_creep->declare_parameters(prm); + diffusion_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); + + std::unique_ptr> dislocation_creep; + dislocation_creep = std::make_unique>(); + dislocation_creep->initialize_simulator (simulator_access.get_simulator()); + dislocation_creep->declare_parameters(prm); + dislocation_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); + + std::unique_ptr> peierls_creep; + peierls_creep = std::make_unique>(); + peierls_creep->initialize_simulator (simulator_access.get_simulator()); + peierls_creep->declare_parameters(prm); + peierls_creep->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); + + std::unique_ptr> drucker_prager_power; + drucker_prager_power = std::make_unique>(); + drucker_prager_power->initialize_simulator (simulator_access.get_simulator()); + drucker_prager_power->declare_parameters(prm); + prm.set("Cohesions", "background:1e9|5e8"); + prm.set("Maximum yield stress", "5e10"); + drucker_prager_power->parse_parameters(prm, std::make_unique>(n_phases_for_each_chemical_composition)); + + // The creep components are arranged in series with each other. + // This package of components is then arranged in parallel with + // a strain rate limiter with a constant viscosity lim_visc. + // The whole system is then arranged in series with a viscosity limiter with + // viscosity max_visc. + // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc) + double min_visc = prm.get_double("Minimum viscosity"); + double max_visc = prm.get_double("Maximum viscosity"); + double lim_visc = (min_visc*max_visc)/(max_visc - min_visc); + + // Assign values to the variables which will be passed to compute_viscosity + // The test involves pure shear calculations at variable pressure and temperature + double temperature; + double pressure; + const double grain_size = 1.e-3; + SymmetricTensor<2,dim> strain_rate; + strain_rate[0][0] = -1e-11; + strain_rate[0][1] = 0.; + strain_rate[1][1] = 1e-11; + strain_rate[2][0] = 0.; + strain_rate[2][1] = 0.; + strain_rate[2][2] = 0.; + + std::cout << "temperature (K) phase transition progress eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl; + + // Loop through strain rates, tracking whether there is a discrepancy in + // the decomposed strain rates. + bool error = false; + double viscosity; + double total_strain_rate; + double creep_strain_rate; + double creep_stress; + double diff_stress; + double disl_stress; + double prls_stress; + double drpr_stress; + std::vector partial_strain_rates(5, 0.); + + for (unsigned int i=0; i <= 2; i++) + { + pressure = 1.e9 + i*2.e9; + std::cout << "pressure: " << pressure / 1.e9 << " GPa" << std::endl; + for (unsigned int j = 0; j <= 10; j++) + { + temperature = 1000. + j*100.; + + // Compute the phase function values + // The depth and gravity are set to zero because they are unused + // when phase functions are calculated by pressure. + // The phase index is set to invalid_unsigned_int, because it is only used internally + // in phase_average_equation_of_state_outputs to loop over all existing phases + MaterialUtilities::PhaseFunctionInputs phase_inputs(temperature, + pressure, + 0., 0., + numbers::invalid_unsigned_int); + + // Compute value of phase functions + for (unsigned int j=0; j < phase_function.n_phase_transitions(); ++j) + { + phase_inputs.phase_transition_index = j; + phase_function_values[j] = phase_function.compute_value(phase_inputs); + } + + Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + + // Compute the viscosity + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates, phase_function_values, n_phase_transitions_per_composition); + total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.); + + // The creep strain rate is calculated by subtracting the strain rate + // of the max viscosity dashpot from the total strain rate + // The creep stress is then calculated by subtracting the stress running + // through the strain rate limiter from the total stress + creep_strain_rate = total_strain_rate - partial_strain_rates[4]; + creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate); + + // Print the output + std::cout << temperature << ' ' << phase_function_values[0] << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate; + for (unsigned int i=0; i < partial_strain_rates.size(); ++i) + { + std::cout << ' ' << partial_strain_rates[i]/total_strain_rate; + } + std::cout << std::endl; + + // The following lines test that each individual creep mechanism + // experiences the same creep stress + + // Each creep mechanism should experience the same stress + diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition); + disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition); + prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition); + if (partial_strain_rates[3] > 0.) + { + drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion, + p.angle_internal_friction, + pressure, + partial_strain_rates[3], + p.max_yield_stress); + } + else + { + drpr_stress = creep_stress; + } + + if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6) + || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6)) + { + error = true; + std::cout << " creep stress: " << creep_stress; + std::cout << " diffusion stress: " << diff_stress; + std::cout << " dislocation stress: " << disl_stress; + std::cout << " peierls stress: " << prls_stress; + std::cout << " drucker prager stress: " << drpr_stress << std::endl; + } + } + } + + if (error) + { + std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl; + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + } + + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } + + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } + + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &) -{ - AssertThrow(false,dealii::ExcInternalError()); -} - -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); -} - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/coordinate_transformation.cc b/tests/coordinate_transformation.cc index 8a5e117df98..2bae9995adb 100644 --- a/tests/coordinate_transformation.cc +++ b/tests/coordinate_transformation.cc @@ -23,71 +23,74 @@ #include -using namespace aspect::Utilities; +namespace aspect +{ + using namespace aspect::Utilities; // Check various conversions between cartesian and spherical coordinates -template -void check_point(T point1, T point2) -{ - std::cout << std::endl << "Point 1: "; - for (unsigned int i = 0; i < dim; ++i) - std::cout << point1[i] << ' '; + template + void check_point(T point1, T point2) + { + std::cout << std::endl << "Point 1: "; + for (unsigned int i = 0; i < dim; ++i) + std::cout << point1[i] << ' '; - std::cout << std::endl << "Point 2: "; - for (unsigned int i = 0; i < dim; ++i) - std::cout << point2[i] << ' '; + std::cout << std::endl << "Point 2: "; + for (unsigned int i = 0; i < dim; ++i) + std::cout << point2[i] << ' '; - std::cout << std::endl; -} + std::cout << std::endl; + } -int f() -{ - const dealii::Point<2> origin2(0,0); - const dealii::Point<3> origin3(0,0,0); + int f() + { + const dealii::Point<2> origin2(0,0); + const dealii::Point<3> origin3(0,0,0); - const std::array sorigin2 = {{0,0}}; - const std::array sorigin3 = {{0,0,0}}; + const std::array sorigin2 = {{0,0}}; + const std::array sorigin3 = {{0,0,0}}; - const dealii::Point<2> one2(1,1); - const dealii::Point<3> one3(1,1,1); + const dealii::Point<2> one2(1,1); + const dealii::Point<3> one3(1,1,1); - const std::array sone2 = {{std::sqrt(2),numbers::PI/4}}; - const std::array sone3 = {{std::sqrt(3),numbers::PI/4,std::acos(1/std::sqrt(3))}}; + const std::array sone2 = {{std::sqrt(2),numbers::PI/4}}; + const std::array sone3 = {{std::sqrt(3),numbers::PI/4,std::acos(1/std::sqrt(3))}}; - const dealii::Point<3> x(1,0,0); - const dealii::Point<3> y(0,1,0); - const dealii::Point<3> z(0,0,1); + const dealii::Point<3> x(1,0,0); + const dealii::Point<3> y(0,1,0); + const dealii::Point<3> z(0,0,1); - const std::array sx = {{1,0,numbers::PI/2}}; - const std::array sy = {{1,numbers::PI/2,numbers::PI/2}}; - const std::array sz = {{1,0,0}}; + const std::array sx = {{1,0,numbers::PI/2}}; + const std::array sy = {{1,numbers::PI/2,numbers::PI/2}}; + const std::array sz = {{1,0,0}}; - check_point,2>(Coordinates::cartesian_to_spherical_coordinates(origin2),sorigin2); - check_point,3>(Coordinates::cartesian_to_spherical_coordinates(origin3),sorigin3); - check_point,2>(origin2, Coordinates::spherical_to_cartesian_coordinates<2>(sorigin2)); - check_point,3>(origin3, Coordinates::spherical_to_cartesian_coordinates<3>(sorigin3)); + check_point,2>(Coordinates::cartesian_to_spherical_coordinates(origin2),sorigin2); + check_point,3>(Coordinates::cartesian_to_spherical_coordinates(origin3),sorigin3); + check_point,2>(origin2, Coordinates::spherical_to_cartesian_coordinates<2>(sorigin2)); + check_point,3>(origin3, Coordinates::spherical_to_cartesian_coordinates<3>(sorigin3)); - check_point,2>(Coordinates::cartesian_to_spherical_coordinates(one2),sone2); - check_point,3>(Coordinates::cartesian_to_spherical_coordinates(one3),sone3); - check_point,2>(one2, Coordinates::spherical_to_cartesian_coordinates<2>(sone2)); - check_point,3>(one3, Coordinates::spherical_to_cartesian_coordinates<3>(sone3)); + check_point,2>(Coordinates::cartesian_to_spherical_coordinates(one2),sone2); + check_point,3>(Coordinates::cartesian_to_spherical_coordinates(one3),sone3); + check_point,2>(one2, Coordinates::spherical_to_cartesian_coordinates<2>(sone2)); + check_point,3>(one3, Coordinates::spherical_to_cartesian_coordinates<3>(sone3)); - check_point,3>(x, Coordinates::spherical_to_cartesian_coordinates<3>(sx)); - check_point,3>(y, Coordinates::spherical_to_cartesian_coordinates<3>(sy)); - check_point,3>(z, Coordinates::spherical_to_cartesian_coordinates<3>(sz)); + check_point,3>(x, Coordinates::spherical_to_cartesian_coordinates<3>(sx)); + check_point,3>(y, Coordinates::spherical_to_cartesian_coordinates<3>(sy)); + check_point,3>(z, Coordinates::spherical_to_cartesian_coordinates<3>(sz)); - check_point,3>(Coordinates::cartesian_to_spherical_coordinates(x),sx); - check_point,3>(Coordinates::cartesian_to_spherical_coordinates(y),sy); - check_point,3>(Coordinates::cartesian_to_spherical_coordinates(z),sz); + check_point,3>(Coordinates::cartesian_to_spherical_coordinates(x),sx); + check_point,3>(Coordinates::cartesian_to_spherical_coordinates(y),sy); + check_point,3>(Coordinates::cartesian_to_spherical_coordinates(z),sz); - const dealii::Point<3> dateline(0,-1,0); - const std::array sdateline = {{1,3*numbers::PI/2,numbers::PI/2}}; + const dealii::Point<3> dateline(0,-1,0); + const std::array sdateline = {{1,3*numbers::PI/2,numbers::PI/2}}; - check_point,3>(dateline, Coordinates::spherical_to_cartesian_coordinates<3>(sdateline)); - check_point,3>(Coordinates::cartesian_to_spherical_coordinates(dateline),sdateline); + check_point,3>(dateline, Coordinates::spherical_to_cartesian_coordinates<3>(sdateline)); + check_point,3>(Coordinates::cartesian_to_spherical_coordinates(dateline),sdateline); - return 42; -} + return 42; + } -int i = f(); + int i = f(); +} diff --git a/tests/drucker_prager_derivatives_2d.cc b/tests/drucker_prager_derivatives_2d.cc index 61404f38416..de7647c1a0e 100644 --- a/tests/drucker_prager_derivatives_2d.cc +++ b/tests/drucker_prager_derivatives_2d.cc @@ -33,214 +33,217 @@ #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - - std::cout << std::endl << "Testing DruckerPrager derivatives against finite difference derivatives " << std::endl; - - using namespace aspect::MaterialModel; - - // first set all material model values - MaterialModelInputs in_base(5,3); - in_base.composition[0][0] = 0; - in_base.composition[0][1] = 0; - in_base.composition[0][2] = 0; - in_base.composition[1][0] = 0.75; - in_base.composition[1][1] = 0.15; - in_base.composition[1][2] = 0.10; - in_base.composition[2][0] = 0; - in_base.composition[2][1] = 0.2; - in_base.composition[2][2] = 0.4; - in_base.composition[3][0] = 0; - in_base.composition[3][1] = 0.2; - in_base.composition[3][2] = 0.4; - in_base.composition[4][0] = 1; - in_base.composition[4][1] = 0; - in_base.composition[4][2] = 0; - - in_base.temperature[0] = 293; - in_base.temperature[1] = 1600; - in_base.temperature[2] = 2000; - in_base.temperature[3] = 2100; - in_base.temperature[4] = 2200; - - in_base.pressure[0] = 1e9; - in_base.pressure[1] = 5e9; - in_base.pressure[2] = 2e10; - in_base.pressure[3] = 2e11; - in_base.pressure[4] = 2e12; - - /** - * We can't take to small strain-rates, because then the difference in the - * visocisty will be too small for the double accuracy which stores - * the visocity solutions and the finite difference solution. - */ - in_base.strain_rate[0] = SymmetricTensor<2,dim>(); - in_base.strain_rate[0][0][0] = 1e-12; - in_base.strain_rate[0][0][1] = 1e-12; - in_base.strain_rate[0][1][1] = 1e-11; - - in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[1][0][0] = -1.71266e-13; - in_base.strain_rate[1][0][1] = -5.82647e-12; - in_base.strain_rate[1][1][1] = 4.21668e-14; - - in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[2][1][1] = 1e-13; - in_base.strain_rate[2][0][1] = 1e-11; - in_base.strain_rate[2][0][0] = -1e-12; - - in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[3][1][1] = 4.9e-21; - in_base.strain_rate[3][0][1] = 4.9e-21; - in_base.strain_rate[3][0][0] = 4.9e-21; - - in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[4][1][1] = 1e-11; - in_base.strain_rate[4][0][1] = 1e-11; - in_base.strain_rate[4][0][0] = 1e-11; - - // initialize some variables we will need later. - const double finite_difference_accuracy = 1e-7; - const double finite_difference_factor = 1+finite_difference_accuracy; - - MaterialModelInputs in_dviscositydstrainrate(in_base); - - MaterialModelOutputs out_base(5,3); - MaterialModelOutputs out_dviscositydpressure(5,3); - MaterialModelOutputs out_dviscositydstrainrate(5,3); - - // initialize the material we want to test. - aspect::ParameterHandler prm; - - const aspect::MaterialModel::DruckerPrager const_material_model = dynamic_cast &>(simulator_access.get_material_model()); - aspect::MaterialModel::DruckerPrager material_model = const_cast &>(const_material_model); - - material_model.declare_parameters(prm); - - prm.enter_subsection("Material model"); + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) { - prm.enter_subsection ("Drucker Prager"); + + std::cout << std::endl << "Testing DruckerPrager derivatives against finite difference derivatives " << std::endl; + + using namespace aspect::MaterialModel; + + // first set all material model values + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 2200; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 2e12; + + /** + * We can't take to small strain-rates, because then the difference in the + * visocisty will be too small for the double accuracy which stores + * the visocity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 4.9e-21; + in_base.strain_rate[3][0][1] = 4.9e-21; + in_base.strain_rate[3][0][0] = 4.9e-21; + + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = 1e-11; + + // initialize some variables we will need later. + const double finite_difference_accuracy = 1e-7; + const double finite_difference_factor = 1+finite_difference_accuracy; + + MaterialModelInputs in_dviscositydstrainrate(in_base); + + MaterialModelOutputs out_base(5,3); + MaterialModelOutputs out_dviscositydpressure(5,3); + MaterialModelOutputs out_dviscositydstrainrate(5,3); + + // initialize the material we want to test. + aspect::ParameterHandler prm; + + const aspect::MaterialModel::DruckerPrager const_material_model = dynamic_cast &>(simulator_access.get_material_model()); + aspect::MaterialModel::DruckerPrager material_model = const_cast &>(const_material_model); + + material_model.declare_parameters(prm); + + prm.enter_subsection("Material model"); { - prm.enter_subsection ("Viscosity"); + prm.enter_subsection ("Drucker Prager"); { - prm.set ("Reference strain rate", "1e-20"); - prm.set ("Angle of internal friction", "30"); + prm.enter_subsection ("Viscosity"); + { + prm.set ("Reference strain rate", "1e-20"); + prm.set ("Angle of internal friction", "30"); + } + prm.leave_subsection(); } prm.leave_subsection(); } prm.leave_subsection(); - } - prm.leave_subsection(); - - const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); - out_base.additional_outputs.push_back(std::make_unique> (5)); + const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); - simulator_access.get_material_model().evaluate(in_base, out_base); + out_base.additional_outputs.push_back(std::make_unique> (5)); - // set up additional output for the derivatives - MaterialModelDerivatives *derivatives; - derivatives = out_base.template get_additional_output>(); - double temp; + simulator_access.get_material_model().evaluate(in_base, out_base); - // have a bool so we know whether the test has succeed or not. - bool Error = false; - - // test the pressure derivative. - MaterialModelInputs in_dviscositydpressure(in_base); - in_dviscositydpressure.pressure[0] *= finite_difference_factor; - in_dviscositydpressure.pressure[1] *= finite_difference_factor; - in_dviscositydpressure.pressure[2] *= finite_difference_factor; - in_dviscositydpressure.pressure[3] *= finite_difference_factor; - in_dviscositydpressure.pressure[4] *= finite_difference_factor; - - simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); - if (in_base.pressure[i] != 0) - { - temp /= (in_base.pressure[i] * finite_difference_accuracy); - } - std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) - { - std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; - Error = true; - } + // set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.template get_additional_output>(); + double temp; - } + // have a bool so we know whether the test has succeed or not. + bool Error = false; - // test the strain-rate derivative. - for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) - { - const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + // test the pressure derivative. + MaterialModelInputs in_dviscositydpressure(in_base); + in_dviscositydpressure.pressure[0] *= finite_difference_factor; + in_dviscositydpressure.pressure[1] *= finite_difference_factor; + in_dviscositydpressure.pressure[2] *= finite_difference_factor; + in_dviscositydpressure.pressure[3] *= finite_difference_factor; + in_dviscositydpressure.pressure[4] *= finite_difference_factor; - for (unsigned int i = 0; i < 5; i++) - { - in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] - + std::fabs(in_base.strain_rate[i][strain_rate_indices]) - * finite_difference_accuracy - * aspect::Utilities::nth_basis_for_symmetric_tensors(component); - } + simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); + if (in_base.pressure[i] != 0) + { + temp /= (in_base.pressure[i] * finite_difference_accuracy); + } + std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) + { + std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; + Error = true; + } - simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + } - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; - } - std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } + // test the strain-rate derivative. + for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) + { + const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + + for (unsigned int i = 0; i < 5; i++) + { + in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] + + std::fabs(in_base.strain_rate[i][strain_rate_indices]) + * finite_difference_accuracy + * aspect::Utilities::nth_basis_for_symmetric_tensors(component); + } + + + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; + } + std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } - } + } - } + if (Error) + { + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } - if (Error) - { - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + } -} + template <> + void f(const aspect::SimulatorAccess<3> &, + aspect::Assemblers::Manager<3> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } -template <> -void f(const aspect::SimulatorAccess<3> &, - aspect::Assemblers::Manager<3> &) -{ - AssertThrow(false,dealii::ExcInternalError()); -} + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/drucker_prager_derivatives_3d.cc b/tests/drucker_prager_derivatives_3d.cc index 2e8c3281178..6ec15235044 100644 --- a/tests/drucker_prager_derivatives_3d.cc +++ b/tests/drucker_prager_derivatives_3d.cc @@ -33,229 +33,232 @@ #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &) +namespace aspect { - - std::cout << std::endl << "Testing DruckerPrager derivatives against finite difference derivatives " << std::endl; - - using namespace aspect::MaterialModel; - - // first set all material model values - MaterialModelInputs in_base(5,3); - in_base.composition[0][0] = 0; - in_base.composition[0][1] = 0; - in_base.composition[0][2] = 0; - in_base.composition[1][0] = 0.75; - in_base.composition[1][1] = 0.15; - in_base.composition[1][2] = 0.10; - in_base.composition[2][0] = 0; - in_base.composition[2][1] = 0.2; - in_base.composition[2][2] = 0.4; - in_base.composition[3][0] = 0; - in_base.composition[3][1] = 0.2; - in_base.composition[3][2] = 0.4; - in_base.composition[4][0] = 1; - in_base.composition[4][1] = 0; - in_base.composition[4][2] = 0; - - in_base.temperature[0] = 293; - in_base.temperature[1] = 1600; - in_base.temperature[2] = 2000; - in_base.temperature[3] = 2100; - in_base.temperature[4] = 2200; - - in_base.pressure[0] = 1e9; - in_base.pressure[1] = 5e9; - in_base.pressure[2] = 2e10; - in_base.pressure[3] = 2e11; - in_base.pressure[4] = 2e12; - - /** - * We can't take to small strain-rates, because then the difference in the - * visocisty will be too small for the double accuracy which stores - * the visocity solutions and the finite difference solution. - */ - in_base.strain_rate[0] = SymmetricTensor<2,dim>(); - in_base.strain_rate[0][0][0] = 1e-12; - in_base.strain_rate[0][0][1] = 1e-12; - in_base.strain_rate[0][1][1] = 1e-11; - in_base.strain_rate[0][2][0] = 1e-12; - in_base.strain_rate[0][2][1] = 1e-12; - in_base.strain_rate[0][2][2] = 1e-11; - - in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[1][0][0] = -1.71266e-13; - in_base.strain_rate[1][0][1] = -5.82647e-12; - in_base.strain_rate[1][1][1] = 4.21668e-14; - in_base.strain_rate[1][2][0] = -5.42647e-12; - in_base.strain_rate[1][2][1] = -5.22647e-12; - in_base.strain_rate[1][2][2] = 4.21668e-14; - - in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[2][1][1] = 1e-13; - in_base.strain_rate[2][0][1] = 1e-11; - in_base.strain_rate[2][0][0] = -1e-12; - in_base.strain_rate[2][2][0] = 1e-11; - in_base.strain_rate[2][2][1] = 1e-11; - in_base.strain_rate[2][2][2] = -1e-12; - - in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[3][1][1] = 4.9e-21; - in_base.strain_rate[3][0][1] = 4.9e-21; - in_base.strain_rate[3][0][0] = 4.9e-21; - in_base.strain_rate[3][2][0] = 4.9e-21; - in_base.strain_rate[3][2][1] = 4.9e-21; - in_base.strain_rate[3][2][2] = 4.9e-21; - - in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[4][1][1] = 1e-11; - in_base.strain_rate[4][0][1] = 1e-11; - in_base.strain_rate[4][0][0] = 1e-11; - in_base.strain_rate[4][2][0] = 1e-11; - in_base.strain_rate[4][2][1] = 1e-11; - in_base.strain_rate[4][2][2] = 1e-11; - - // initialize some variables we will need later. - const double finite_difference_accuracy = 1e-7; - const double finite_difference_factor = 1+finite_difference_accuracy; - - MaterialModelInputs in_dviscositydstrainrate(in_base); - - MaterialModelOutputs out_base(5,3); - MaterialModelOutputs out_dviscositydpressure(5,3); - MaterialModelOutputs out_dviscositydstrainrate(5,3); - - // initialize the material we want to test. - aspect::ParameterHandler prm; - - const aspect::MaterialModel::DruckerPrager const_material_model = dynamic_cast &>(simulator_access.get_material_model()); - aspect::MaterialModel::DruckerPrager material_model = const_cast &>(const_material_model); - - material_model.declare_parameters(prm); - - prm.enter_subsection("Material model"); + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &) { - prm.enter_subsection ("Drucker Prager"); + + std::cout << std::endl << "Testing DruckerPrager derivatives against finite difference derivatives " << std::endl; + + using namespace aspect::MaterialModel; + + // first set all material model values + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 2200; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 2e12; + + /** + * We can't take to small strain-rates, because then the difference in the + * visocisty will be too small for the double accuracy which stores + * the visocity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + in_base.strain_rate[0][2][0] = 1e-12; + in_base.strain_rate[0][2][1] = 1e-12; + in_base.strain_rate[0][2][2] = 1e-11; + + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + in_base.strain_rate[1][2][0] = -5.42647e-12; + in_base.strain_rate[1][2][1] = -5.22647e-12; + in_base.strain_rate[1][2][2] = 4.21668e-14; + + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + in_base.strain_rate[2][2][0] = 1e-11; + in_base.strain_rate[2][2][1] = 1e-11; + in_base.strain_rate[2][2][2] = -1e-12; + + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 4.9e-21; + in_base.strain_rate[3][0][1] = 4.9e-21; + in_base.strain_rate[3][0][0] = 4.9e-21; + in_base.strain_rate[3][2][0] = 4.9e-21; + in_base.strain_rate[3][2][1] = 4.9e-21; + in_base.strain_rate[3][2][2] = 4.9e-21; + + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = 1e-11; + in_base.strain_rate[4][2][0] = 1e-11; + in_base.strain_rate[4][2][1] = 1e-11; + in_base.strain_rate[4][2][2] = 1e-11; + + // initialize some variables we will need later. + const double finite_difference_accuracy = 1e-7; + const double finite_difference_factor = 1+finite_difference_accuracy; + + MaterialModelInputs in_dviscositydstrainrate(in_base); + + MaterialModelOutputs out_base(5,3); + MaterialModelOutputs out_dviscositydpressure(5,3); + MaterialModelOutputs out_dviscositydstrainrate(5,3); + + // initialize the material we want to test. + aspect::ParameterHandler prm; + + const aspect::MaterialModel::DruckerPrager const_material_model = dynamic_cast &>(simulator_access.get_material_model()); + aspect::MaterialModel::DruckerPrager material_model = const_cast &>(const_material_model); + + material_model.declare_parameters(prm); + + prm.enter_subsection("Material model"); { - prm.enter_subsection ("Viscosity"); + prm.enter_subsection ("Drucker Prager"); { - prm.set ("Reference strain rate", "1e-20"); - prm.set ("Angle of internal friction", "30"); + prm.enter_subsection ("Viscosity"); + { + prm.set ("Reference strain rate", "1e-20"); + prm.set ("Angle of internal friction", "30"); + } + prm.leave_subsection(); } prm.leave_subsection(); } prm.leave_subsection(); - } - prm.leave_subsection(); - - const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); - out_base.additional_outputs.push_back(std::make_unique> (5)); + const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); - simulator_access.get_material_model().evaluate(in_base, out_base); + out_base.additional_outputs.push_back(std::make_unique> (5)); - // set up additional output for the derivatives - MaterialModelDerivatives *derivatives; - derivatives = out_base.template get_additional_output>(); - double temp; + simulator_access.get_material_model().evaluate(in_base, out_base); - // have a bool so we know whether the test has succeed or not. - bool Error = false; - - // test the pressure derivative. - MaterialModelInputs in_dviscositydpressure(in_base); - in_dviscositydpressure.pressure[0] *= finite_difference_factor; - in_dviscositydpressure.pressure[1] *= finite_difference_factor; - in_dviscositydpressure.pressure[2] *= finite_difference_factor; - in_dviscositydpressure.pressure[3] *= finite_difference_factor; - in_dviscositydpressure.pressure[4] *= finite_difference_factor; - - simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); - if (in_base.pressure[i] != 0) - { - temp /= (in_base.pressure[i] * finite_difference_accuracy); - } - std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) - { - std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; - Error = true; - } + // set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.template get_additional_output>(); + double temp; - } + // have a bool so we know whether the test has succeed or not. + bool Error = false; - // test the strain-rate derivative. - for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) - { - const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + // test the pressure derivative. + MaterialModelInputs in_dviscositydpressure(in_base); + in_dviscositydpressure.pressure[0] *= finite_difference_factor; + in_dviscositydpressure.pressure[1] *= finite_difference_factor; + in_dviscositydpressure.pressure[2] *= finite_difference_factor; + in_dviscositydpressure.pressure[3] *= finite_difference_factor; + in_dviscositydpressure.pressure[4] *= finite_difference_factor; - for (unsigned int i = 0; i < 5; i++) - { - in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] - + std::fabs(in_base.strain_rate[i][strain_rate_indices]) - * finite_difference_accuracy - * aspect::Utilities::nth_basis_for_symmetric_tensors(component); - } + simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); + if (in_base.pressure[i] != 0) + { + temp /= (in_base.pressure[i] * finite_difference_accuracy); + } + std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) + { + std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; + Error = true; + } - simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + } - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; - } - std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } + // test the strain-rate derivative. + for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) + { + const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + + for (unsigned int i = 0; i < 5; i++) + { + in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] + + std::fabs(in_base.strain_rate[i][strain_rate_indices]) + * finite_difference_accuracy + * aspect::Utilities::nth_basis_for_symmetric_tensors(component); + } + + + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; + } + std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } - } + } - } + if (Error) + { + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } - if (Error) - { - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + } -} + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &) + { + AssertThrow(false,dealii::ExcInternalError()); + } -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &) -{ - AssertThrow(false,dealii::ExcInternalError()); -} + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2)); + } -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2)); + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/own_gravity.cc b/tests/own_gravity.cc index e4ef8057cb6..aac5dd918cb 100644 --- a/tests/own_gravity.cc +++ b/tests/own_gravity.cc @@ -22,20 +22,23 @@ #include #include -template -class MyGravity : - public aspect::GravityModel::Interface +namespace aspect { - public: - virtual Tensor<1,dim> gravity_vector (const Point &position) const - { - Tensor<1,dim> ret; - ret[0] = position[1]; - ret[1] = 42.0; - return ret; - } -}; + template + class MyGravity : + public aspect::GravityModel::Interface + { + public: + virtual Tensor<1,dim> gravity_vector (const Point &position) const + { + Tensor<1,dim> ret; + ret[0] = position[1]; + ret[1] = 42.0; + return ret; + } + }; // explicit instantiation -ASPECT_REGISTER_GRAVITY_MODEL(MyGravity, "my gravity", "no description") + ASPECT_REGISTER_GRAVITY_MODEL(MyGravity, "my gravity", "no description") +} diff --git a/tests/prm_distance_polygon.cc b/tests/prm_distance_polygon.cc index 50d48fda0f8..85bcdf1dde5 100644 --- a/tests/prm_distance_polygon.cc +++ b/tests/prm_distance_polygon.cc @@ -24,54 +24,57 @@ #include -int f() +namespace aspect { - using namespace aspect::Utilities; - - const int dim=3; - - // A square polygon - std::vector> polygon(4); - polygon[0] = Point<2>(0.0,0.0); - polygon[1] = Point<2>(1.0,0.0); - polygon[2] = Point<2>(1.0,1.0); - polygon[3] = Point<2>(0.0,1.0); - - // A concave polygon - std::vector> concave_polygon(5); - concave_polygon[0] = Point<2>(0.0,0.0); - concave_polygon[1] = Point<2>(1.0,0.0); - concave_polygon[2] = Point<2>(1.0,1.0); - concave_polygon[3] = Point<2>(0.5,0.5); - concave_polygon[4] = Point<2>(0.0,1.0); - - // A selfcrossing polygon - std::vector> crossing_polygon(4); - crossing_polygon[0] = Point<2>(0.0,0.0); - crossing_polygon[1] = Point<2>(1.0,0.0); - crossing_polygon[2] = Point<2>(1.0,-1.0); - crossing_polygon[3] = Point<2>(0.0,1.0); - - - // Some points inside and outside the polygon - Point<2> points[] = {Point<2>(0.5,-1), Point<2>(0.5,0.5), Point<2>(0.001,0.2), Point<2>(2.0,2.0), Point<2>(0.25,0.70)}; - - std::cout << "Testing distance to polygon function with the following parameters: (polygon 1) " - << polygon[0] << ", " << polygon[1] << ", " << polygon[2] << ", " << polygon[3] << ", " - << "(polygon 2) " << concave_polygon[0] << ", " << concave_polygon[1] << ", " << concave_polygon[2] << ", " << concave_polygon[3] << ", " << concave_polygon[4] << ", " - << "(polygon 3) " << crossing_polygon[0] << ", " << crossing_polygon[1] << ", " << crossing_polygon[2] << ", " << crossing_polygon[3] - << ", (points) " - << points[0] << ", " << points[1] << ", " << points[2] << ", " << points[3] << ", " << points[4] << std::endl; - - for (unsigned int i = 0; i < 5; i++) - { - std::cout << "Minimal distance of point " << points[i] << " to polygon 1 = " << signed_distance_to_polygon(polygon,points[i]) << std::endl; - std::cout << "Minimal distance of point " << points[i] << " to polygon 2 = " << signed_distance_to_polygon(concave_polygon,points[i]) << std::endl; - std::cout << "Minimal distance of point " << points[i] << " to polygon 3 = " << signed_distance_to_polygon(crossing_polygon,points[i]) << std::endl; - } - - exit(0); - return 42; -} + int f() + { + using namespace aspect::Utilities; + + const int dim=3; + + // A square polygon + std::vector> polygon(4); + polygon[0] = Point<2>(0.0,0.0); + polygon[1] = Point<2>(1.0,0.0); + polygon[2] = Point<2>(1.0,1.0); + polygon[3] = Point<2>(0.0,1.0); + + // A concave polygon + std::vector> concave_polygon(5); + concave_polygon[0] = Point<2>(0.0,0.0); + concave_polygon[1] = Point<2>(1.0,0.0); + concave_polygon[2] = Point<2>(1.0,1.0); + concave_polygon[3] = Point<2>(0.5,0.5); + concave_polygon[4] = Point<2>(0.0,1.0); + + // A selfcrossing polygon + std::vector> crossing_polygon(4); + crossing_polygon[0] = Point<2>(0.0,0.0); + crossing_polygon[1] = Point<2>(1.0,0.0); + crossing_polygon[2] = Point<2>(1.0,-1.0); + crossing_polygon[3] = Point<2>(0.0,1.0); + + + // Some points inside and outside the polygon + Point<2> points[] = {Point<2>(0.5,-1), Point<2>(0.5,0.5), Point<2>(0.001,0.2), Point<2>(2.0,2.0), Point<2>(0.25,0.70)}; + + std::cout << "Testing distance to polygon function with the following parameters: (polygon 1) " + << polygon[0] << ", " << polygon[1] << ", " << polygon[2] << ", " << polygon[3] << ", " + << "(polygon 2) " << concave_polygon[0] << ", " << concave_polygon[1] << ", " << concave_polygon[2] << ", " << concave_polygon[3] << ", " << concave_polygon[4] << ", " + << "(polygon 3) " << crossing_polygon[0] << ", " << crossing_polygon[1] << ", " << crossing_polygon[2] << ", " << crossing_polygon[3] + << ", (points) " + << points[0] << ", " << points[1] << ", " << points[2] << ", " << points[3] << ", " << points[4] << std::endl; + + for (unsigned int i = 0; i < 5; i++) + { + std::cout << "Minimal distance of point " << points[i] << " to polygon 1 = " << signed_distance_to_polygon(polygon,points[i]) << std::endl; + std::cout << "Minimal distance of point " << points[i] << " to polygon 2 = " << signed_distance_to_polygon(concave_polygon,points[i]) << std::endl; + std::cout << "Minimal distance of point " << points[i] << " to polygon 3 = " << signed_distance_to_polygon(crossing_polygon,points[i]) << std::endl; + } + + exit(0); + return 42; + } // run this function by initializing a global variable by it -int i = f(); + int i = f(); +} diff --git a/tests/simple_nonlinear.cc b/tests/simple_nonlinear.cc index 023770ed4c5..c16f554b64f 100644 --- a/tests/simple_nonlinear.cc +++ b/tests/simple_nonlinear.cc @@ -30,224 +30,228 @@ #include "../benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/simple_nonlinear.cc" -template -int f(double parameter) +namespace aspect { + template + int f(double parameter) + { - std::cout << std::endl << "Test for p = " << parameter << " with dimension " << dim << std::endl; - - using namespace aspect::MaterialModel; - - // first set all material model values - MaterialModelInputs in_base(5,3); - in_base.composition[0][0] = 0; - in_base.composition[0][1] = 0; - in_base.composition[0][2] = 0; - in_base.composition[1][0] = 0.75; - in_base.composition[1][1] = 0.15; - in_base.composition[1][2] = 0.10; - in_base.composition[2][0] = 0; - in_base.composition[2][1] = 0.2; - in_base.composition[2][2] = 0.4; - in_base.composition[3][0] = 0; - in_base.composition[3][1] = 0.2; - in_base.composition[3][2] = 0.4; - in_base.composition[4][0] = 1; - in_base.composition[4][1] = 0; - in_base.composition[4][2] = 0; - - in_base.temperature[0] = 293; - in_base.temperature[1] = 1600; - in_base.temperature[2] = 2000; - in_base.temperature[3] = 2100; - in_base.temperature[4] = 2200; - - in_base.pressure[0] = 1e9; - in_base.pressure[1] = 5e9; - in_base.pressure[2] = 2e10; - in_base.pressure[3] = 2e11; - in_base.pressure[4] = 2e12; - - /** - * We can't take to small strain-rates, because then the difference in the - * viscosity will be too small for the double accuracy which stores - * the viscosity solutions and the finite difference solution. - */ - in_base.strain_rate[0] = SymmetricTensor<2,dim>(); - in_base.strain_rate[0][0][0] = 1e-12; - in_base.strain_rate[0][0][1] = 1e-12; - in_base.strain_rate[0][1][1] = 1e-11; - if (dim == 3) - { - in_base.strain_rate[0][2][0] = 1e-12; - in_base.strain_rate[0][2][1] = 1e-12; - in_base.strain_rate[0][2][2] = 1e-11; - } - - in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[1][0][0] = -1.71266e-13; - in_base.strain_rate[1][0][1] = -5.82647e-12; - in_base.strain_rate[1][1][1] = 4.21668e-14; - if (dim == 3) - { - in_base.strain_rate[1][2][0] = -5.42647e-12; - in_base.strain_rate[1][2][1] = -5.22647e-12; - in_base.strain_rate[1][2][2] = 4.21668e-14; - } - in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[2][1][1] = 1e-13; - in_base.strain_rate[2][0][1] = 1e-11; - in_base.strain_rate[2][0][0] = -1e-12; - if (dim == 3) + std::cout << std::endl << "Test for p = " << parameter << " with dimension " << dim << std::endl; + + using namespace aspect::MaterialModel; + + // first set all material model values + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 2200; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 2e12; + + /** + * We can't take to small strain-rates, because then the difference in the + * viscosity will be too small for the double accuracy which stores + * the viscosity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + if (dim == 3) + { + in_base.strain_rate[0][2][0] = 1e-12; + in_base.strain_rate[0][2][1] = 1e-12; + in_base.strain_rate[0][2][2] = 1e-11; + } + + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + if (dim == 3) + { + in_base.strain_rate[1][2][0] = -5.42647e-12; + in_base.strain_rate[1][2][1] = -5.22647e-12; + in_base.strain_rate[1][2][2] = 4.21668e-14; + } + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + if (dim == 3) + { + in_base.strain_rate[2][2][0] = 1e-11; + in_base.strain_rate[2][2][1] = 1e-11; + in_base.strain_rate[2][2][2] = -1e-12; + } + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 4.9e-21; + in_base.strain_rate[3][0][1] = 4.9e-21; + in_base.strain_rate[3][0][0] = 4.9e-21; + if (dim == 3) + { + in_base.strain_rate[3][2][0] = 4.9e-21; + in_base.strain_rate[3][2][1] = 4.9e-21; + in_base.strain_rate[3][2][2] = 4.9e-21; + } + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = 1e-11; + if (dim == 3) + { + in_base.strain_rate[4][2][0] = 1e-11; + in_base.strain_rate[4][2][1] = 1e-11; + in_base.strain_rate[4][2][2] = 1e-11; + } + + // initialize some variables we will need later. + double finite_difference_accuracy = 1e-7; + double finite_difference_factor = 1+finite_difference_accuracy; + + + MaterialModelInputs in_dviscositydstrainrate(in_base); + + MaterialModelOutputs out_base(5,3); + MaterialModelOutputs out_dviscositydstrainrate(5,3); + + if (out_base.template get_additional_output>() != nullptr) + throw "error"; + + out_base.additional_outputs.push_back(std::make_unique> (5)); + + // initialize the material we want to test. + SimpleNonlinear mat; + ParameterHandler prm; + mat.declare_parameters(prm); + + prm.enter_subsection("Compositional fields"); { - in_base.strain_rate[2][2][0] = 1e-11; - in_base.strain_rate[2][2][1] = 1e-11; - in_base.strain_rate[2][2][2] = -1e-12; + prm.set("Number of fields","3"); } - in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[3][1][1] = 4.9e-21; - in_base.strain_rate[3][0][1] = 4.9e-21; - in_base.strain_rate[3][0][0] = 4.9e-21; - if (dim == 3) - { - in_base.strain_rate[3][2][0] = 4.9e-21; - in_base.strain_rate[3][2][1] = 4.9e-21; - in_base.strain_rate[3][2][2] = 4.9e-21; - } - in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[4][1][1] = 1e-11; - in_base.strain_rate[4][0][1] = 1e-11; - in_base.strain_rate[4][0][0] = 1e-11; - if (dim == 3) + prm.leave_subsection(); + prm.enter_subsection("Material model"); { - in_base.strain_rate[4][2][0] = 1e-11; - in_base.strain_rate[4][2][1] = 1e-11; - in_base.strain_rate[4][2][2] = 1e-11; + prm.enter_subsection ("Simple nonlinear"); + { + prm.set ("Viscosity prefactor", "1e-37,1e-36,1e-35,5e-36"); + prm.set ("Viscosity averaging p", std::to_string(parameter)); + prm.set ("Minimum strain rate", 1.4e-20); + } + prm.leave_subsection(); } + prm.leave_subsection(); - // initialize some variables we will need later. - double finite_difference_accuracy = 1e-7; - double finite_difference_factor = 1+finite_difference_accuracy; - - - MaterialModelInputs in_dviscositydstrainrate(in_base); - - MaterialModelOutputs out_base(5,3); - MaterialModelOutputs out_dviscositydstrainrate(5,3); - - if (out_base.template get_additional_output>() != nullptr) - throw "error"; - - out_base.additional_outputs.push_back(std::make_unique> (5)); - - // initialize the material we want to test. - SimpleNonlinear mat; - ParameterHandler prm; - mat.declare_parameters(prm); - - prm.enter_subsection("Compositional fields"); - { - prm.set("Number of fields","3"); + mat.parse_parameters(prm); + + mat.evaluate(in_base, out_base); + + // set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.template get_additional_output>(); + double temp; + + // have a bool so we know whether the test has succeed or not. + bool Error = false; + + // this material is not pressure dependent, so we do not test it. + + // test the strain-rate derivative. + for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) + { + const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + + for (unsigned int i = 0; i < 5; i++) + { + in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] + + std::fabs(in_base.strain_rate[i][strain_rate_indices]) + * finite_difference_accuracy + * aspect::Utilities::nth_basis_for_symmetric_tensors(component); + } + + + mat.evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; + } + std::cout << "strain-rate: component = " << component << ", point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } + + } + + if (Error) + { + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + + return 42; } - prm.leave_subsection(); - prm.enter_subsection("Material model"); + + int exit_function() { - prm.enter_subsection ("Simple nonlinear"); - { - prm.set ("Viscosity prefactor", "1e-37,1e-36,1e-35,5e-36"); - prm.set ("Viscosity averaging p", std::to_string(parameter)); - prm.set ("Minimum strain rate", 1.4e-20); - } - prm.leave_subsection(); + exit(0); + return 42; } - prm.leave_subsection(); - - mat.parse_parameters(prm); - - mat.evaluate(in_base, out_base); - - // set up additional output for the derivatives - MaterialModelDerivatives *derivatives; - derivatives = out_base.template get_additional_output>(); - double temp; - - // have a bool so we know whether the test has succeed or not. - bool Error = false; - - // this material is not pressure dependent, so we do not test it. - - // test the strain-rate derivative. - for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) - { - const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); - - for (unsigned int i = 0; i < 5; i++) - { - in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] - + std::fabs(in_base.strain_rate[i][strain_rate_indices]) - * finite_difference_accuracy - * aspect::Utilities::nth_basis_for_symmetric_tensors(component); - } - - - mat.evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; - } - std::cout << "strain-rate: component = " << component << ", point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } - - } - - } - - if (Error) - { - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } - - return 42; -} - -int exit_function() -{ - exit(0); - return 42; -} // run this function by initializing a global variable by it // test 2D -int ii2 = f<2>(-1000); // Testing min function -int iz2 = f<2>(-2); // Testing generalized p norm mean with negative p -int ij2 = f<2>(-1.5); // Testing generalized p norm mean with negative, non int p -int ik2 = f<2>(-1); // Testing harmonic mean -int ji2 = f<2>(0); // Testing geometric mean -int jj2 = f<2>(1); // Testing arithmetic mean -int jk2 = f<2>(2); // Testing generalized p norm mean with positive p -int kj2 = f<2>(1000); // Testing max function + int ii2 = f<2>(-1000); // Testing min function + int iz2 = f<2>(-2); // Testing generalized p norm mean with negative p + int ij2 = f<2>(-1.5); // Testing generalized p norm mean with negative, non int p + int ik2 = f<2>(-1); // Testing harmonic mean + int ji2 = f<2>(0); // Testing geometric mean + int jj2 = f<2>(1); // Testing arithmetic mean + int jk2 = f<2>(2); // Testing generalized p norm mean with positive p + int kj2 = f<2>(1000); // Testing max function // test 3D -int ii3 = f<3>(-1000); // Testing min function -int iz3 = f<3>(-2); // Testing generalized p norm mean with negative p -int ij3 = f<3>(-1.5); // Testing generalized p norm mean with negative, non int p -int ik3 = f<3>(-1); // Testing harmonic mean -int ji3 = f<3>(0); // Testing geometric mean -int jj3 = f<3>(1); // Testing arithmetic mean -int jk3 = f<3>(2); // Testing generalized p norm mean with positive p -int kj3 = f<3>(1000); // Testing max function + int ii3 = f<3>(-1000); // Testing min function + int iz3 = f<3>(-2); // Testing generalized p norm mean with negative p + int ij3 = f<3>(-1.5); // Testing generalized p norm mean with negative, non int p + int ik3 = f<3>(-1); // Testing harmonic mean + int ji3 = f<3>(0); // Testing geometric mean + int jj3 = f<3>(1); // Testing arithmetic mean + int jk3 = f<3>(2); // Testing generalized p norm mean with positive p + int kj3 = f<3>(1000); // Testing max function // exit -int kl2 = exit_function(); + int kl2 = exit_function(); + +} diff --git a/tests/spiegelman_material.cc b/tests/spiegelman_material.cc index 5650b030bd8..c4ddb9280c6 100644 --- a/tests/spiegelman_material.cc +++ b/tests/spiegelman_material.cc @@ -20,284 +20,288 @@ #include "../tests/spiegelman_fail_test.cc" -int f(double parameter) +namespace aspect { - - std::cout << std::endl << "Test for p = " << parameter << std::endl; - - const int dim=2; - using namespace aspect::MaterialModel; - MaterialModelInputs in_base(5,3); - in_base.composition[0][0] = 0; - in_base.composition[0][1] = 0; - in_base.composition[0][2] = 0; - in_base.composition[1][0] = 0.75; - in_base.composition[1][1] = 0.15; - in_base.composition[1][2] = 0.10; - in_base.composition[2][0] = 0; - in_base.composition[2][1] = 0.2; - in_base.composition[2][2] = 0.4; - in_base.composition[3][0] = 0; - in_base.composition[3][1] = 0.2; - in_base.composition[3][2] = 0.4; - in_base.composition[4][0] = 1; - in_base.composition[4][1] = 0; - in_base.composition[4][2] = 0; - - in_base.pressure[0] = 1e9; - in_base.pressure[1] = 5e9; - in_base.pressure[2] = 2e10; - in_base.pressure[3] = 2e11; - in_base.pressure[4] = 2e12; - - /** - * We can't take to small strain-rates, because then the difference in the - * viscosity will be too small for the double accuracy which stores - * the viscosity solutions and the finite difference solution. - */ - in_base.strain_rate[0] = SymmetricTensor<2,dim>(); - in_base.strain_rate[0][0][0] = 1e-12; - in_base.strain_rate[0][0][1] = 1e-12; - in_base.strain_rate[0][1][1] = 1e-11; - in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[1][0][0] = -1.71266e-13; - in_base.strain_rate[1][0][1] = -5.82647e-12; - in_base.strain_rate[1][1][1] = 4.21668e-14; - in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[2][1][1] = 1e-13; - in_base.strain_rate[2][0][1] = 1e-11; - in_base.strain_rate[2][0][0] = -1e-12; - in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[3][1][1] = 1e-22; - in_base.strain_rate[3][0][1] = 1e-22; - in_base.strain_rate[3][0][0] = -1e-22; - /** - * We can get it working with these values: - * - * in_base.strain_rate[3][1][1] = 9e-16; - * in_base.strain_rate[3][0][1] = 9e-16; - * in_base.strain_rate[3][0][0] = -9e-16; - * - * With very low strain-rates the differences between the finite difference and - * analytical solution grow. We interpent this as that the finite difference - * method we use becomes unaccurate for very low strain-rates. - */ - in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[4][1][1] = 1e-11; - in_base.strain_rate[4][0][1] = 1e-11; - in_base.strain_rate[4][0][0] = 1e-11; - - in_base.temperature[0] = 293; - in_base.temperature[1] = 1600; - in_base.temperature[2] = 2000; - in_base.temperature[3] = 2100; - in_base.temperature[4] = 2200; - - SymmetricTensor<2,dim> zerozero = SymmetricTensor<2,dim>(); - SymmetricTensor<2,dim> onezero = SymmetricTensor<2,dim>(); - SymmetricTensor<2,dim> oneone = SymmetricTensor<2,dim>(); - - zerozero[0][0] = 1; - onezero[1][0] = 0.5; // because symmetry doubles this entry - oneone[1][1] = 1; - - double finite_difference_accuracy = 1e-7; - double finite_difference_factor = 1+finite_difference_accuracy; - - bool Error = false; - - MaterialModelInputs in_dviscositydpressure(in_base); - in_dviscositydpressure.pressure[0] *= finite_difference_factor; - in_dviscositydpressure.pressure[1] *= finite_difference_factor; - in_dviscositydpressure.pressure[2] *= finite_difference_factor; - in_dviscositydpressure.pressure[3] *= finite_difference_factor; - in_dviscositydpressure.pressure[4] *= finite_difference_factor; - - MaterialModelInputs in_dviscositydstrainrate_zerozero(in_base); - MaterialModelInputs in_dviscositydstrainrate_onezero(in_base); - MaterialModelInputs in_dviscositydstrainrate_oneone(in_base); - - in_dviscositydstrainrate_zerozero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[0][0][0]) * finite_difference_accuracy * zerozero; - in_dviscositydstrainrate_zerozero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[1][0][0]) * finite_difference_accuracy * zerozero; - in_dviscositydstrainrate_zerozero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[2][0][0]) * finite_difference_accuracy * zerozero; - in_dviscositydstrainrate_zerozero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[3][0][0]) * finite_difference_accuracy * zerozero; - in_dviscositydstrainrate_zerozero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[4][0][0]) * finite_difference_accuracy * zerozero; - in_dviscositydstrainrate_onezero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[0][1][0]) * finite_difference_accuracy * onezero; - in_dviscositydstrainrate_onezero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[1][1][0]) * finite_difference_accuracy * onezero; - in_dviscositydstrainrate_onezero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[2][1][0]) * finite_difference_accuracy * onezero; - in_dviscositydstrainrate_onezero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[3][1][0]) * finite_difference_accuracy * onezero; - in_dviscositydstrainrate_onezero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[4][1][0]) * finite_difference_accuracy * onezero; - in_dviscositydstrainrate_oneone.strain_rate[0] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[0][1][1]) * finite_difference_accuracy * oneone; - in_dviscositydstrainrate_oneone.strain_rate[1] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[1][1][1]) * finite_difference_accuracy * oneone; - in_dviscositydstrainrate_oneone.strain_rate[2] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[2][1][1]) * finite_difference_accuracy * oneone; - in_dviscositydstrainrate_oneone.strain_rate[3] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[3][1][1]) * finite_difference_accuracy * oneone; - in_dviscositydstrainrate_oneone.strain_rate[4] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[4][1][1]) * finite_difference_accuracy * oneone; - - MaterialModelInputs in_dviscositydtemperature(in_base); - in_dviscositydtemperature.temperature[0] *= 1.0000000001; - in_dviscositydtemperature.temperature[1] *= 1.0000000001; - in_dviscositydtemperature.temperature[2] *= 1.0000000001; - in_dviscositydtemperature.temperature[3] *= 1.0000000001; - in_dviscositydtemperature.temperature[4] *= 1.0000000001; - - - MaterialModelOutputs out_base(5,3); - - MaterialModelOutputs out_dviscositydpressure(5,3); - MaterialModelOutputs out_dviscositydstrainrate_zerozero(5,3); - MaterialModelOutputs out_dviscositydstrainrate_onezero(5,3); - MaterialModelOutputs out_dviscositydstrainrate_oneone(5,3); - MaterialModelOutputs out_dviscositydtemperature(5,3); - - if (out_base.get_additional_output>() != nullptr) - throw "error"; - - out_base.additional_outputs.push_back(std::make_unique> (5)); - - SpiegelmanMaterial mat; - ParameterHandler prm; - mat.declare_parameters(prm); - - prm.enter_subsection("Compositional fields"); - { - prm.set("Number of fields","3"); - prm.set("List of conductivities","2.25"); - prm.set("List of capacities","1250"); - prm.set("List of reference densities","2700.0"); - prm.set("List of cohesions","1e8,0,1e8,1e8"); - prm.set("List of angles of internal friction","30.0"); - prm.set("List of initial viscosities","1e20"); - prm.set("List of constant viscosities","0,1e21,0,0"); - } - prm.leave_subsection(); - prm.enter_subsection("Material model"); + int f(double parameter) { - prm.enter_subsection ("Spiegelman 2016"); - { - prm.set("Use deviator of strain-rate", "false"); - prm.set("Use analytical derivative", "true"); - prm.set ("Viscosity averaging p", std::to_string(parameter)); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - - mat.parse_parameters(prm); - - mat.evaluate(in_base, out_base); - mat.evaluate(in_dviscositydpressure, out_dviscositydpressure); - mat.evaluate(in_dviscositydstrainrate_zerozero, out_dviscositydstrainrate_zerozero); - mat.evaluate(in_dviscositydstrainrate_onezero, out_dviscositydstrainrate_onezero); - mat.evaluate(in_dviscositydstrainrate_oneone, out_dviscositydstrainrate_oneone); - mat.evaluate(in_dviscositydtemperature, out_dviscositydtemperature); - //set up additional output for the derivatives - MaterialModelDerivatives *derivatives; - derivatives = out_base.get_additional_output>(); - - double temp; - for (unsigned int i = 0; i < 5; i++) + std::cout << std::endl << "Test for p = " << parameter << std::endl; + + const int dim=2; + using namespace aspect::MaterialModel; + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 2e12; + + /** + * We can't take to small strain-rates, because then the difference in the + * viscosity will be too small for the double accuracy which stores + * the viscosity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 1e-22; + in_base.strain_rate[3][0][1] = 1e-22; + in_base.strain_rate[3][0][0] = -1e-22; + /** + * We can get it working with these values: + * + * in_base.strain_rate[3][1][1] = 9e-16; + * in_base.strain_rate[3][0][1] = 9e-16; + * in_base.strain_rate[3][0][0] = -9e-16; + * + * With very low strain-rates the differences between the finite difference and + * analytical solution grow. We interpent this as that the finite difference + * method we use becomes unaccurate for very low strain-rates. + */ + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = 1e-11; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 2200; + + SymmetricTensor<2,dim> zerozero = SymmetricTensor<2,dim>(); + SymmetricTensor<2,dim> onezero = SymmetricTensor<2,dim>(); + SymmetricTensor<2,dim> oneone = SymmetricTensor<2,dim>(); + + zerozero[0][0] = 1; + onezero[1][0] = 0.5; // because symmetry doubles this entry + oneone[1][1] = 1; + + double finite_difference_accuracy = 1e-7; + double finite_difference_factor = 1+finite_difference_accuracy; + + bool Error = false; + + MaterialModelInputs in_dviscositydpressure(in_base); + in_dviscositydpressure.pressure[0] *= finite_difference_factor; + in_dviscositydpressure.pressure[1] *= finite_difference_factor; + in_dviscositydpressure.pressure[2] *= finite_difference_factor; + in_dviscositydpressure.pressure[3] *= finite_difference_factor; + in_dviscositydpressure.pressure[4] *= finite_difference_factor; + + MaterialModelInputs in_dviscositydstrainrate_zerozero(in_base); + MaterialModelInputs in_dviscositydstrainrate_onezero(in_base); + MaterialModelInputs in_dviscositydstrainrate_oneone(in_base); + + in_dviscositydstrainrate_zerozero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[0][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[1][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[2][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[3][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_zerozero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[4][0][0]) * finite_difference_accuracy * zerozero; + in_dviscositydstrainrate_onezero.strain_rate[0] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[0][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[1] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[1][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[2] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[2][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[3] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[3][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_onezero.strain_rate[4] += std::fabs(in_dviscositydstrainrate_onezero.strain_rate[4][1][0]) * finite_difference_accuracy * onezero; + in_dviscositydstrainrate_oneone.strain_rate[0] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[0][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[1] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[1][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[2] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[2][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[3] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[3][1][1]) * finite_difference_accuracy * oneone; + in_dviscositydstrainrate_oneone.strain_rate[4] += std::fabs(in_dviscositydstrainrate_oneone.strain_rate[4][1][1]) * finite_difference_accuracy * oneone; + + MaterialModelInputs in_dviscositydtemperature(in_base); + in_dviscositydtemperature.temperature[0] *= 1.0000000001; + in_dviscositydtemperature.temperature[1] *= 1.0000000001; + in_dviscositydtemperature.temperature[2] *= 1.0000000001; + in_dviscositydtemperature.temperature[3] *= 1.0000000001; + in_dviscositydtemperature.temperature[4] *= 1.0000000001; + + + MaterialModelOutputs out_base(5,3); + + MaterialModelOutputs out_dviscositydpressure(5,3); + MaterialModelOutputs out_dviscositydstrainrate_zerozero(5,3); + MaterialModelOutputs out_dviscositydstrainrate_onezero(5,3); + MaterialModelOutputs out_dviscositydstrainrate_oneone(5,3); + MaterialModelOutputs out_dviscositydtemperature(5,3); + + if (out_base.get_additional_output>() != nullptr) + throw "error"; + + out_base.additional_outputs.push_back(std::make_unique> (5)); + + SpiegelmanMaterial mat; + ParameterHandler prm; + mat.declare_parameters(prm); + + prm.enter_subsection("Compositional fields"); { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); - if (in_base.pressure[i] != 0) - { - temp /= (in_base.pressure[i] * finite_difference_accuracy); - } - - std::cout << "pressure on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) - { - std::cout << "Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; - Error = true; - } - + prm.set("Number of fields","3"); + prm.set("List of conductivities","2.25"); + prm.set("List of capacities","1250"); + prm.set("List of reference densities","2700.0"); + prm.set("List of cohesions","1e8,0,1e8,1e8"); + prm.set("List of angles of internal friction","30.0"); + prm.set("List of initial viscosities","1e20"); + prm.set("List of constant viscosities","0,1e21,0,0"); } - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate_zerozero.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[i][0][0]) * finite_difference_accuracy; - } - std::cout << "zerozero on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][0][0] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][0][0]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][0][0]))) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } - - - - } - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate_onezero.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate_onezero.strain_rate[i][1][0]) * finite_difference_accuracy; - } - std::cout << "onezero on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][0] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][0]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][0])) ) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } - } - - for (unsigned int i = 0; i < 5; i++) + prm.leave_subsection(); + prm.enter_subsection("Material model"); { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate_oneone.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate_oneone.strain_rate[i][1][1]) * finite_difference_accuracy; - } - std::cout << "oneone on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][1] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][1]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][1])) ) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } - + prm.enter_subsection ("Spiegelman 2016"); + { + prm.set("Use deviator of strain-rate", "false"); + prm.set("Use analytical derivative", "true"); + prm.set ("Viscosity averaging p", std::to_string(parameter)); + } + prm.leave_subsection(); } + prm.leave_subsection(); - if (Error) - { - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + mat.parse_parameters(prm); + + mat.evaluate(in_base, out_base); + mat.evaluate(in_dviscositydpressure, out_dviscositydpressure); + mat.evaluate(in_dviscositydstrainrate_zerozero, out_dviscositydstrainrate_zerozero); + mat.evaluate(in_dviscositydstrainrate_onezero, out_dviscositydstrainrate_onezero); + mat.evaluate(in_dviscositydstrainrate_oneone, out_dviscositydstrainrate_oneone); + mat.evaluate(in_dviscositydtemperature, out_dviscositydtemperature); + + //set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.get_additional_output>(); + + double temp; + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); + if (in_base.pressure[i] != 0) + { + temp /= (in_base.pressure[i] * finite_difference_accuracy); + } + + std::cout << "pressure on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) + { + std::cout << "Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; + Error = true; + } + + } + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate_zerozero.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate_zerozero.strain_rate[i][0][0]) * finite_difference_accuracy; + } + std::cout << "zerozero on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][0][0] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][0][0]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][0][0]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + + + } + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate_onezero.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate_onezero.strain_rate[i][1][0]) * finite_difference_accuracy; + } + std::cout << "onezero on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][0] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][0]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][0])) ) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + } + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate_oneone.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate_oneone.strain_rate[i][1][1]) * finite_difference_accuracy; + } + std::cout << "oneone on quadrature point " << i << ": Finite difference = " << temp << ". Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][1][1] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][1][1]) > 1e-3 * 0.5 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][1][1])) ) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } + + if (Error) + { + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } + + return 42; + } - return 42; -} + int exit_function() + { + exit(0); + return 42; + } +// run this function by initializing a global variable by it + int ii = f(-1000); // Testing min function + int iz = f(-2); // Testing generalized p norm mean with negative p + int ij = f(-1.5); // Testing generalized p norm mean with negative, non int p + int ik = f(-1); // Testing harmonic mean + int ji = f(0); // Testing geometric mean + int jj = f(1); // Testing arithmetic mean + int jk = f(2); // Testing generalized p norm mean with positive p + int kj = f(1000); // Testing max function + int kl = exit_function(); -int exit_function() -{ - exit(0); - return 42; } -// run this function by initializing a global variable by it -int ii = f(-1000); // Testing min function -int iz = f(-2); // Testing generalized p norm mean with negative p -int ij = f(-1.5); // Testing generalized p norm mean with negative, non int p -int ik = f(-1); // Testing harmonic mean -int ji = f(0); // Testing geometric mean -int jj = f(1); // Testing arithmetic mean -int jk = f(2); // Testing generalized p norm mean with positive p -int kj = f(1000); // Testing max function -int kl = exit_function(); diff --git a/tests/visco_plastic_derivatives_2d.cc b/tests/visco_plastic_derivatives_2d.cc index 1468d4563ff..f97b46d729c 100644 --- a/tests/visco_plastic_derivatives_2d.cc +++ b/tests/visco_plastic_derivatives_2d.cc @@ -34,240 +34,243 @@ #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &, - std::string averaging_parameter) +namespace aspect { - - std::cout << std::endl << "Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter " << averaging_parameter << std::endl; - - using namespace aspect::MaterialModel; - - // first set all material model values - MaterialModelInputs in_base(5,3); - in_base.composition[0][0] = 0; - in_base.composition[0][1] = 0; - in_base.composition[0][2] = 0; - in_base.composition[1][0] = 0.75; - in_base.composition[1][1] = 0.15; - in_base.composition[1][2] = 0.10; - in_base.composition[2][0] = 0; - in_base.composition[2][1] = 0.2; - in_base.composition[2][2] = 0.4; - in_base.composition[3][0] = 0; - in_base.composition[3][1] = 0.2; - in_base.composition[3][2] = 0.4; - in_base.composition[4][0] = 1; - in_base.composition[4][1] = 0; - in_base.composition[4][2] = 0; - - in_base.temperature[0] = 293; - in_base.temperature[1] = 1600; - in_base.temperature[2] = 2000; - in_base.temperature[3] = 2100; - in_base.temperature[4] = 600; - - in_base.pressure[0] = 1e9; - in_base.pressure[1] = 5e9; - in_base.pressure[2] = 2e10; - in_base.pressure[3] = 2e11; - in_base.pressure[4] = 5e8; - - in_base.position[0] = Point(); - in_base.position[1] = Point(); - in_base.position[2] = Point(); - in_base.position[3] = Point(); - in_base.position[4] = Point(); - - /** - * We can't take too small strain-rates, because then the difference in the - * viscosity will be too small for the double accuracy which stores - * the viscosity solutions and the finite difference solution. - */ - in_base.strain_rate[0] = SymmetricTensor<2,dim>(); - in_base.strain_rate[0][0][0] = 1e-12; - in_base.strain_rate[0][0][1] = 1e-12; - in_base.strain_rate[0][1][1] = 1e-11; - - in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[1][0][0] = -1.71266e-13; - in_base.strain_rate[1][0][1] = -5.82647e-12; - in_base.strain_rate[1][1][1] = 4.21668e-14; - - in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[2][1][1] = 1e-13; - in_base.strain_rate[2][0][1] = 1e-11; - in_base.strain_rate[2][0][0] = -1e-12; - - in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[3][1][1] = 4.9e-21; - in_base.strain_rate[3][0][1] = 4.9e-21; - in_base.strain_rate[3][0][0] = 4.9e-21; - - in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[4][1][1] = 1e-11; - in_base.strain_rate[4][0][1] = 1e-11; - in_base.strain_rate[4][0][0] = -1e-11; - - - // initialize some variables we will need later. - double finite_difference_accuracy = 1e-7; - double finite_difference_factor = 1+finite_difference_accuracy; - - - MaterialModelInputs in_dviscositydstrainrate(in_base); - - MaterialModelOutputs out_base(5,3); - MaterialModelOutputs out_dviscositydpressure(5,3); - MaterialModelOutputs out_dviscositydstrainrate(5,3); - - // initialize the material we want to test. - aspect::ParameterHandler prm; - - const aspect::MaterialModel::ViscoPlastic &const_material_model = dynamic_cast &>(simulator_access.get_material_model()); - aspect::MaterialModel::ViscoPlastic &material_model = const_cast &>(const_material_model); - - material_model.declare_parameters(prm); - - prm.enter_subsection("Material model"); + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &, + std::string averaging_parameter) { - prm.enter_subsection ("Visco Plastic"); - { - prm.set ("Viscosity averaging scheme", averaging_parameter); - prm.set ("Angles of internal friction", "30"); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - - const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); - - out_base.additional_outputs.push_back(std::make_unique> (5)); - - simulator_access.get_material_model().evaluate(in_base, out_base); - // set up additional output for the derivatives - MaterialModelDerivatives *derivatives; - derivatives = out_base.template get_additional_output>(); - double temp; - - // have a bool so we know whether the test has succeed or not. - bool Error = false; - - // test the pressure derivative. - MaterialModelInputs in_dviscositydpressure(in_base); - in_dviscositydpressure.pressure[0] *= finite_difference_factor; - in_dviscositydpressure.pressure[1] *= finite_difference_factor; - in_dviscositydpressure.pressure[2] *= finite_difference_factor; - in_dviscositydpressure.pressure[3] *= finite_difference_factor; - in_dviscositydpressure.pressure[4] *= finite_difference_factor; - - simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); - - for (unsigned int i = 0; i < 5; i++) + std::cout << std::endl << "Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter " << averaging_parameter << std::endl; + + using namespace aspect::MaterialModel; + + // first set all material model values + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 600; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 5e8; + + in_base.position[0] = Point(); + in_base.position[1] = Point(); + in_base.position[2] = Point(); + in_base.position[3] = Point(); + in_base.position[4] = Point(); + + /** + * We can't take too small strain-rates, because then the difference in the + * viscosity will be too small for the double accuracy which stores + * the viscosity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 4.9e-21; + in_base.strain_rate[3][0][1] = 4.9e-21; + in_base.strain_rate[3][0][0] = 4.9e-21; + + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = -1e-11; + + + // initialize some variables we will need later. + double finite_difference_accuracy = 1e-7; + double finite_difference_factor = 1+finite_difference_accuracy; + + + MaterialModelInputs in_dviscositydstrainrate(in_base); + + MaterialModelOutputs out_base(5,3); + MaterialModelOutputs out_dviscositydpressure(5,3); + MaterialModelOutputs out_dviscositydstrainrate(5,3); + + // initialize the material we want to test. + aspect::ParameterHandler prm; + + const aspect::MaterialModel::ViscoPlastic &const_material_model = dynamic_cast &>(simulator_access.get_material_model()); + aspect::MaterialModel::ViscoPlastic &material_model = const_cast &>(const_material_model); + + material_model.declare_parameters(prm); + + prm.enter_subsection("Material model"); { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); - if (in_base.pressure[i] != 0) - { - temp /= (in_base.pressure[i] * finite_difference_accuracy); - } - std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) - { - std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; - Error = true; - } - + prm.enter_subsection ("Visco Plastic"); + { + prm.set ("Viscosity averaging scheme", averaging_parameter); + prm.set ("Angles of internal friction", "30"); + } + prm.leave_subsection(); } + prm.leave_subsection(); - // test the strain-rate derivative. - for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) - { - const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); - - for (unsigned int i = 0; i < 5; i++) - { - // components that are not on the diagonal are multiplied by 0.5, because the symmetric tensor - // is modified by 0.5 in both symmetric directions (xy/yx) simultaneously and we compute the combined - // derivative - in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] - + std::fabs(in_base.strain_rate[i][strain_rate_indices]) - * (component > dim-1 ? 0.5 : 1 ) - * finite_difference_accuracy - * aspect::Utilities::nth_basis_for_symmetric_tensors(component); - } - - - simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; - } - std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } - - } + const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); + + out_base.additional_outputs.push_back(std::make_unique> (5)); + + simulator_access.get_material_model().evaluate(in_base, out_base); + + // set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.template get_additional_output>(); + double temp; + + // have a bool so we know whether the test has succeed or not. + bool Error = false; + + // test the pressure derivative. + MaterialModelInputs in_dviscositydpressure(in_base); + in_dviscositydpressure.pressure[0] *= finite_difference_factor; + in_dviscositydpressure.pressure[1] *= finite_difference_factor; + in_dviscositydpressure.pressure[2] *= finite_difference_factor; + in_dviscositydpressure.pressure[3] *= finite_difference_factor; + in_dviscositydpressure.pressure[4] *= finite_difference_factor; + + simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); + if (in_base.pressure[i] != 0) + { + temp /= (in_base.pressure[i] * finite_difference_accuracy); + } + std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) + { + std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; + Error = true; + } + + } + + // test the strain-rate derivative. + for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) + { + const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + + for (unsigned int i = 0; i < 5; i++) + { + // components that are not on the diagonal are multiplied by 0.5, because the symmetric tensor + // is modified by 0.5 in both symmetric directions (xy/yx) simultaneously and we compute the combined + // derivative + in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] + + std::fabs(in_base.strain_rate[i][strain_rate_indices]) + * (component > dim-1 ? 0.5 : 1 ) + * finite_difference_accuracy + * aspect::Utilities::nth_basis_for_symmetric_tensors(component); + } + + + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; + } + std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } + + } + + if (Error) + { + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } - } - - if (Error) - { - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + } -} + template <> + void f(const aspect::SimulatorAccess<3> &, + aspect::Assemblers::Manager<3> &, + std::string ) + { + AssertThrow(false,dealii::ExcInternalError()); + } -template <> -void f(const aspect::SimulatorAccess<3> &, - aspect::Assemblers::Manager<3> &, - std::string ) -{ - AssertThrow(false,dealii::ExcInternalError()); -} + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "harmonic")); + + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "geometric")); + + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "arithmetic")); + + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "maximum composition")); + } -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "harmonic")); - - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "geometric")); - - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "arithmetic")); - - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "maximum composition")); + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>) diff --git a/tests/visco_plastic_derivatives_3d.cc b/tests/visco_plastic_derivatives_3d.cc index dcf32291bcd..930a6a4ec51 100644 --- a/tests/visco_plastic_derivatives_3d.cc +++ b/tests/visco_plastic_derivatives_3d.cc @@ -34,255 +34,258 @@ #include -template -void f(const aspect::SimulatorAccess &simulator_access, - aspect::Assemblers::Manager &, - std::string averaging_parameter) +namespace aspect { - - std::cout << std::endl << "Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter " << averaging_parameter << std::endl; - - using namespace aspect::MaterialModel; - - // first set all material model values - MaterialModelInputs in_base(5,3); - in_base.composition[0][0] = 0; - in_base.composition[0][1] = 0; - in_base.composition[0][2] = 0; - in_base.composition[1][0] = 0.75; - in_base.composition[1][1] = 0.15; - in_base.composition[1][2] = 0.10; - in_base.composition[2][0] = 0; - in_base.composition[2][1] = 0.2; - in_base.composition[2][2] = 0.4; - in_base.composition[3][0] = 0; - in_base.composition[3][1] = 0.2; - in_base.composition[3][2] = 0.4; - in_base.composition[4][0] = 1; - in_base.composition[4][1] = 0; - in_base.composition[4][2] = 0; - - in_base.temperature[0] = 293; - in_base.temperature[1] = 1600; - in_base.temperature[2] = 2000; - in_base.temperature[3] = 2100; - in_base.temperature[4] = 600; - - in_base.pressure[0] = 1e9; - in_base.pressure[1] = 5e9; - in_base.pressure[2] = 2e10; - in_base.pressure[3] = 2e11; - in_base.pressure[4] = 5e8; - - in_base.position[0] = Point(); - in_base.position[1] = Point(); - in_base.position[2] = Point(); - in_base.position[3] = Point(); - in_base.position[4] = Point(); - - /** - * We can't take too small strain-rates, because then the difference in the - * viscosity will be too small for the double accuracy which stores - * the viscosity solutions and the finite difference solution. - */ - in_base.strain_rate[0] = SymmetricTensor<2,dim>(); - in_base.strain_rate[0][0][0] = 1e-12; - in_base.strain_rate[0][0][1] = 1e-12; - in_base.strain_rate[0][1][1] = 1e-11; - in_base.strain_rate[0][2][0] = 1e-12; - in_base.strain_rate[0][2][1] = 1e-12; - in_base.strain_rate[0][2][2] = 1e-11; - - in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[1][0][0] = -1.71266e-13; - in_base.strain_rate[1][0][1] = -5.82647e-12; - in_base.strain_rate[1][1][1] = 4.21668e-14; - in_base.strain_rate[1][2][0] = -5.42647e-12; - in_base.strain_rate[1][2][1] = -5.22647e-12; - in_base.strain_rate[1][2][2] = 4.21668e-14; - - in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[2][1][1] = 1e-13; - in_base.strain_rate[2][0][1] = 1e-11; - in_base.strain_rate[2][0][0] = -1e-12; - in_base.strain_rate[2][2][0] = 1e-11; - in_base.strain_rate[2][2][1] = 1e-11; - in_base.strain_rate[2][2][2] = -1e-12; - - in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[3][1][1] = 4.9e-21; - in_base.strain_rate[3][0][1] = 4.9e-21; - in_base.strain_rate[3][0][0] = 4.9e-21; - in_base.strain_rate[3][2][0] = 4.9e-21; - in_base.strain_rate[3][2][1] = 4.9e-21; - in_base.strain_rate[3][2][2] = 4.9e-21; - - in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); - in_base.strain_rate[4][1][1] = 1e-11; - in_base.strain_rate[4][0][1] = 1e-11; - in_base.strain_rate[4][0][0] = -5e-12; - in_base.strain_rate[4][2][0] = 1e-11; - in_base.strain_rate[4][2][1] = 1e-11; - in_base.strain_rate[4][2][2] = -5e-12; - - - // initialize some variables we will need later. - double finite_difference_accuracy = 1e-7; - double finite_difference_factor = 1+finite_difference_accuracy; - - - MaterialModelInputs in_dviscositydstrainrate(in_base); - - MaterialModelOutputs out_base(5,3); - MaterialModelOutputs out_dviscositydpressure(5,3); - MaterialModelOutputs out_dviscositydstrainrate(5,3); - - // initialize the material we want to test. - aspect::ParameterHandler prm; - - const aspect::MaterialModel::ViscoPlastic &const_material_model = dynamic_cast &>(simulator_access.get_material_model()); - aspect::MaterialModel::ViscoPlastic &material_model = const_cast &>(const_material_model); - - material_model.declare_parameters(prm); - - prm.enter_subsection("Material model"); + template + void f(const aspect::SimulatorAccess &simulator_access, + aspect::Assemblers::Manager &, + std::string averaging_parameter) { - prm.enter_subsection ("Visco Plastic"); - { - prm.set ("Viscosity averaging scheme", averaging_parameter); - prm.set ("Angles of internal friction", "30"); - } - prm.leave_subsection(); - } - prm.leave_subsection(); - - const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); - - out_base.additional_outputs.push_back(std::make_unique> (5)); - - simulator_access.get_material_model().evaluate(in_base, out_base); - // set up additional output for the derivatives - MaterialModelDerivatives *derivatives; - derivatives = out_base.template get_additional_output>(); - double temp; - - // have a bool so we know whether the test has succeed or not. - bool Error = false; - - // test the pressure derivative. - MaterialModelInputs in_dviscositydpressure(in_base); - in_dviscositydpressure.pressure[0] *= finite_difference_factor; - in_dviscositydpressure.pressure[1] *= finite_difference_factor; - in_dviscositydpressure.pressure[2] *= finite_difference_factor; - in_dviscositydpressure.pressure[3] *= finite_difference_factor; - in_dviscositydpressure.pressure[4] *= finite_difference_factor; - - simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); - - for (unsigned int i = 0; i < 5; i++) + std::cout << std::endl << "Testing ViscoPlastic derivatives against analytical derivatives for averaging parameter " << averaging_parameter << std::endl; + + using namespace aspect::MaterialModel; + + // first set all material model values + MaterialModelInputs in_base(5,3); + in_base.composition[0][0] = 0; + in_base.composition[0][1] = 0; + in_base.composition[0][2] = 0; + in_base.composition[1][0] = 0.75; + in_base.composition[1][1] = 0.15; + in_base.composition[1][2] = 0.10; + in_base.composition[2][0] = 0; + in_base.composition[2][1] = 0.2; + in_base.composition[2][2] = 0.4; + in_base.composition[3][0] = 0; + in_base.composition[3][1] = 0.2; + in_base.composition[3][2] = 0.4; + in_base.composition[4][0] = 1; + in_base.composition[4][1] = 0; + in_base.composition[4][2] = 0; + + in_base.temperature[0] = 293; + in_base.temperature[1] = 1600; + in_base.temperature[2] = 2000; + in_base.temperature[3] = 2100; + in_base.temperature[4] = 600; + + in_base.pressure[0] = 1e9; + in_base.pressure[1] = 5e9; + in_base.pressure[2] = 2e10; + in_base.pressure[3] = 2e11; + in_base.pressure[4] = 5e8; + + in_base.position[0] = Point(); + in_base.position[1] = Point(); + in_base.position[2] = Point(); + in_base.position[3] = Point(); + in_base.position[4] = Point(); + + /** + * We can't take too small strain-rates, because then the difference in the + * viscosity will be too small for the double accuracy which stores + * the viscosity solutions and the finite difference solution. + */ + in_base.strain_rate[0] = SymmetricTensor<2,dim>(); + in_base.strain_rate[0][0][0] = 1e-12; + in_base.strain_rate[0][0][1] = 1e-12; + in_base.strain_rate[0][1][1] = 1e-11; + in_base.strain_rate[0][2][0] = 1e-12; + in_base.strain_rate[0][2][1] = 1e-12; + in_base.strain_rate[0][2][2] = 1e-11; + + in_base.strain_rate[1] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[1][0][0] = -1.71266e-13; + in_base.strain_rate[1][0][1] = -5.82647e-12; + in_base.strain_rate[1][1][1] = 4.21668e-14; + in_base.strain_rate[1][2][0] = -5.42647e-12; + in_base.strain_rate[1][2][1] = -5.22647e-12; + in_base.strain_rate[1][2][2] = 4.21668e-14; + + in_base.strain_rate[2] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[2][1][1] = 1e-13; + in_base.strain_rate[2][0][1] = 1e-11; + in_base.strain_rate[2][0][0] = -1e-12; + in_base.strain_rate[2][2][0] = 1e-11; + in_base.strain_rate[2][2][1] = 1e-11; + in_base.strain_rate[2][2][2] = -1e-12; + + in_base.strain_rate[3] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[3][1][1] = 4.9e-21; + in_base.strain_rate[3][0][1] = 4.9e-21; + in_base.strain_rate[3][0][0] = 4.9e-21; + in_base.strain_rate[3][2][0] = 4.9e-21; + in_base.strain_rate[3][2][1] = 4.9e-21; + in_base.strain_rate[3][2][2] = 4.9e-21; + + in_base.strain_rate[4] = SymmetricTensor<2,dim>(in_base.strain_rate[0]); + in_base.strain_rate[4][1][1] = 1e-11; + in_base.strain_rate[4][0][1] = 1e-11; + in_base.strain_rate[4][0][0] = -5e-12; + in_base.strain_rate[4][2][0] = 1e-11; + in_base.strain_rate[4][2][1] = 1e-11; + in_base.strain_rate[4][2][2] = -5e-12; + + + // initialize some variables we will need later. + double finite_difference_accuracy = 1e-7; + double finite_difference_factor = 1+finite_difference_accuracy; + + + MaterialModelInputs in_dviscositydstrainrate(in_base); + + MaterialModelOutputs out_base(5,3); + MaterialModelOutputs out_dviscositydpressure(5,3); + MaterialModelOutputs out_dviscositydstrainrate(5,3); + + // initialize the material we want to test. + aspect::ParameterHandler prm; + + const aspect::MaterialModel::ViscoPlastic &const_material_model = dynamic_cast &>(simulator_access.get_material_model()); + aspect::MaterialModel::ViscoPlastic &material_model = const_cast &>(const_material_model); + + material_model.declare_parameters(prm); + + prm.enter_subsection("Material model"); { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); - if (in_base.pressure[i] != 0) - { - temp /= (in_base.pressure[i] * finite_difference_accuracy); - } - std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) - { - std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; - Error = true; - } - + prm.enter_subsection ("Visco Plastic"); + { + prm.set ("Viscosity averaging scheme", averaging_parameter); + prm.set ("Angles of internal friction", "30"); + } + prm.leave_subsection(); } + prm.leave_subsection(); - // test the strain-rate derivative. - for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) - { - const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); - - for (unsigned int i = 0; i < 5; i++) - { - // components that are not on the diagonal are multiplied by 0.5, because the symmetric tensor - // is modified by 0.5 in both symmetric directions (xy/yx) simultaneously and we compute the combined - // derivative - in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] - + std::fabs(in_base.strain_rate[i][strain_rate_indices]) - * (component > dim-1 ? 0.5 : 1 ) - * finite_difference_accuracy - * aspect::Utilities::nth_basis_for_symmetric_tensors(component); - } - - - simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); - - for (unsigned int i = 0; i < 5; i++) - { - // prevent division by zero. If it is zero, the test has passed, because or - // the finite difference and the analytical result match perfectly, or (more - // likely) the material model in independent of this variable. - temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; - if (temp != 0) - { - temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; - } - std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; - if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) - { - std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; - Error = true; - } - - } + const_cast &>(simulator_access.get_material_model()).parse_parameters(prm); + + out_base.additional_outputs.push_back(std::make_unique> (5)); + + simulator_access.get_material_model().evaluate(in_base, out_base); + + // set up additional output for the derivatives + MaterialModelDerivatives *derivatives; + derivatives = out_base.template get_additional_output>(); + double temp; + + // have a bool so we know whether the test has succeed or not. + bool Error = false; + + // test the pressure derivative. + MaterialModelInputs in_dviscositydpressure(in_base); + in_dviscositydpressure.pressure[0] *= finite_difference_factor; + in_dviscositydpressure.pressure[1] *= finite_difference_factor; + in_dviscositydpressure.pressure[2] *= finite_difference_factor; + in_dviscositydpressure.pressure[3] *= finite_difference_factor; + in_dviscositydpressure.pressure[4] *= finite_difference_factor; + + simulator_access.get_material_model().evaluate(in_dviscositydpressure, out_dviscositydpressure); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = (out_dviscositydpressure.viscosities[i] - out_base.viscosities[i]); + if (in_base.pressure[i] != 0) + { + temp /= (in_base.pressure[i] * finite_difference_accuracy); + } + std::cout << "pressure: point = " << i << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_pressure[i] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_pressure[i]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_pressure[i]))) + { + std::cout << " Error: The derivative of the viscosity to the pressure is too different from the analytical value." << std::endl; + Error = true; + } + + } + + // test the strain-rate derivative. + for (unsigned int component = 0; component < SymmetricTensor<2,dim>::n_independent_components; ++component) + { + const TableIndices<2> strain_rate_indices = SymmetricTensor<2,dim>::unrolled_to_component_indices (component); + + for (unsigned int i = 0; i < 5; i++) + { + // components that are not on the diagonal are multiplied by 0.5, because the symmetric tensor + // is modified by 0.5 in both symmetric directions (xy/yx) simultaneously and we compute the combined + // derivative + in_dviscositydstrainrate.strain_rate[i] = in_base.strain_rate[i] + + std::fabs(in_base.strain_rate[i][strain_rate_indices]) + * (component > dim-1 ? 0.5 : 1 ) + * finite_difference_accuracy + * aspect::Utilities::nth_basis_for_symmetric_tensors(component); + } + + + simulator_access.get_material_model().evaluate(in_dviscositydstrainrate, out_dviscositydstrainrate); + + for (unsigned int i = 0; i < 5; i++) + { + // prevent division by zero. If it is zero, the test has passed, because or + // the finite difference and the analytical result match perfectly, or (more + // likely) the material model in independent of this variable. + temp = out_dviscositydstrainrate.viscosities[i] - out_base.viscosities[i]; + if (temp != 0) + { + temp /= std::fabs(in_dviscositydstrainrate.strain_rate[i][strain_rate_indices]) * finite_difference_accuracy; + } + std::cout << "strain-rate: point = " << i << ", component = " << component << ", Finite difference = " << temp << ", Analytical derivative = " << derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices] << std::endl; + if (std::fabs(temp - derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]) > 1e-3 * (std::fabs(temp) + std::fabs(derivatives->viscosity_derivative_wrt_strain_rate[i][strain_rate_indices]))) + { + std::cout << " Error: The derivative of the viscosity to the strain rate is too different from the analytical value." << std::endl; + Error = true; + } + + } + + } + + if (Error) + { + std::cout << "Some parts of the test were not successful." << std::endl; + } + else + { + std::cout << "OK" << std::endl; + } - } - - if (Error) - { - std::cout << "Some parts of the test were not successful." << std::endl; - } - else - { - std::cout << "OK" << std::endl; - } + } -} + template <> + void f(const aspect::SimulatorAccess<2> &, + aspect::Assemblers::Manager<2> &, + std::string ) + { + AssertThrow(false,dealii::ExcInternalError()); + } -template <> -void f(const aspect::SimulatorAccess<2> &, - aspect::Assemblers::Manager<2> &, - std::string ) -{ - AssertThrow(false,dealii::ExcInternalError()); -} + template + void signal_connector (aspect::SimulatorSignals &signals) + { + std::cout << "* Connecting signals" << std::endl; + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "harmonic")); + + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "geometric")); + + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "arithmetic")); + + signals.set_assemblers.connect (std::bind(&f, + std::placeholders::_1, + std::placeholders::_2, + "maximum composition")); + } -template -void signal_connector (aspect::SimulatorSignals &signals) -{ - std::cout << "* Connecting signals" << std::endl; - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "harmonic")); - - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "geometric")); - - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "arithmetic")); - - signals.set_assemblers.connect (std::bind(&f, - std::placeholders::_1, - std::placeholders::_2, - "maximum composition")); + ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) } - -ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, - signal_connector<3>)