From c6366cef5bf15d685e92f657d72a21bc1e8c6d9d Mon Sep 17 00:00:00 2001 From: Jonathan Robey Date: Wed, 3 Aug 2016 14:35:02 -0700 Subject: [PATCH 1/5] Remove duplication of constraint computation in initialization --- source/simulator/initial_conditions.cc | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index 549dbfe0488..ebbc709455f 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -51,6 +51,10 @@ namespace aspect // we need to track whether we need to normalize the totality of fields bool normalize_composition = false; + // we need the current constraints to be correct for + // the current time + compute_current_constraints (); + #if DEAL_II_VERSION_GTE(8,5,0) initial_solution.reinit(system_rhs, false); @@ -70,9 +74,6 @@ namespace aspect // then apply constraints and copy the // result into vectors with ghost elements. to do so, - // we need the current constraints to be correct for - // the current time - compute_current_constraints (); current_constraints.distribute(initial_solution); // copy temperature/composition block only @@ -163,9 +164,6 @@ namespace aspect // then apply constraints and copy the // result into vectors with ghost elements. to do so, - // we need the current constraints to be correct for - // the current time - compute_current_constraints (); current_constraints.distribute(initial_solution); // copy temperature/composition block only @@ -283,9 +281,6 @@ namespace aspect // then apply constraints and copy the // result into vectors with ghost elements. to do so, - // we need the current constraints to be correct for - // the current time - compute_current_constraints (); current_constraints.distribute(initial_solution); // copy temperature/composition block only From 3cc745a617eb9286f4c55dd35b8d5c7615b9e8a0 Mon Sep 17 00:00:00 2001 From: Jonathan Robey Date: Wed, 3 Aug 2016 15:27:57 -0700 Subject: [PATCH 2/5] Add test for composition normalization --- tests/normalize_initial_composition.prm | 61 +++++++++++++++++++ .../screen-output | 24 ++++++++ 2 files changed, 85 insertions(+) create mode 100644 tests/normalize_initial_composition.prm create mode 100644 tests/normalize_initial_composition/screen-output diff --git a/tests/normalize_initial_composition.prm b/tests/normalize_initial_composition.prm new file mode 100644 index 00000000000..d728a039b51 --- /dev/null +++ b/tests/normalize_initial_composition.prm @@ -0,0 +1,61 @@ +# This test ensures that the Simpler material model works +# with compositional fields enabled + +set Dimension = 2 +set End time = 0 + +subsection Boundary temperature model + set Model name = initial temperature +end + +subsection Compositional fields + set Number of fields = 2 + set List of normalized fields = 0,1 +end + +subsection Compositional initial conditions + subsection Function + set Function expression = 1;1 + end +end + +subsection Geometry model + set Model name = box + subsection Box + set X extent = 1 + set X repetitions = 2 + set Y extent = 1 + set Y repetitions = 2 + end +end + +subsection Gravity model + set Model name = vertical + subsection Vertical + set Magnitude = 1 + end +end + +subsection Initial conditions + set Model name = function + subsection Function + set Function expression = 1 + end +end + +subsection Material model + set Model name = simple +end + +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 +end + +subsection Model settings + set Tangential velocity boundary indicators = left, right, top, bottom +end + +subsection Postprocess + set List of postprocessors = composition statistics +end diff --git a/tests/normalize_initial_composition/screen-output b/tests/normalize_initial_composition/screen-output new file mode 100644 index 00000000000..5c849b3798c --- /dev/null +++ b/tests/normalize_initial_composition/screen-output @@ -0,0 +1,24 @@ +----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +Number of active cells: 4 (on 1 levels) +Number of degrees of freedom: 134 (50+9+25+25+25) + +Sum of compositional fields is not one, fields will be normalized +*** Timestep 0: t=0 years + Solving temperature system... 0 iterations. + Solving C_1 system ... 0 iterations. + Solving C_2 system ... 0 iterations. + Rebuilding Stokes preconditioner... + Solving Stokes system... 2+0 iterations. + + Postprocessing: + Compositions min/max/mass: 0.5/0.5/0.5 // 0.5/0.5/0.5 + +Termination requested by criterion: end time + + ++---------------------------------------------+------------+------------+ ++---------------------------------+-----------+------------+------------+ ++---------------------------------+-----------+------------+------------+ + From 209b56bc6f6b22c027d01b56d994e4799b9f6b67 Mon Sep 17 00:00:00 2001 From: Jonathan Robey Date: Wed, 3 Aug 2016 16:08:58 -0700 Subject: [PATCH 3/5] Remove constraint/copy duplication --- source/simulator/initial_conditions.cc | 64 +++++++++----------------- 1 file changed, 22 insertions(+), 42 deletions(-) diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index ebbc709455f..f97f2319018 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -51,12 +51,9 @@ namespace aspect // we need to track whether we need to normalize the totality of fields bool normalize_composition = false; - // we need the current constraints to be correct for - // the current time - compute_current_constraints (); + initial_solution.reinit(system_rhs, false); #if DEAL_II_VERSION_GTE(8,5,0) - initial_solution.reinit(system_rhs, false); // Interpolate initial temperature VectorFunctionFromScalarFunctionObject @@ -72,21 +69,10 @@ namespace aspect initial_solution, introspection.component_masks.temperature); - // then apply constraints and copy the - // result into vectors with ghost elements. to do so, - current_constraints.distribute(initial_solution); - - // copy temperature/composition block only - const unsigned int t_blockidx = introspection.block_indices.temperature; - solution.block(t_blockidx) = initial_solution.block(t_blockidx); - old_solution.block(t_blockidx) = initial_solution.block(t_blockidx); - old_old_solution.block(t_blockidx) = initial_solution.block(t_blockidx); - // Interpolate initial compositions for (unsigned int n=0; n composition_function (std_cxx11::bind(&CompositionalInitialConditions::Interface::initial_composition, @@ -105,7 +91,7 @@ namespace aspect const unsigned int base_element = advf.base_element(introspection); - // get the temperature/composition support points + // get the composition support points const std::vector > support_points = finite_element.base_element(base_element).get_unit_support_points(); Assert (support_points.size() != 0, @@ -129,7 +115,7 @@ namespace aspect { // if it is specified in the parameter file that the sum of all compositional fields // must not exceed one, this should be checked - if (parameters.normalized_fields.size()>0 && n == 1) + if (parameters.normalized_fields.size()>0 && n == 0) { double sum = 0; for (unsigned int m=0; m Date: Wed, 3 Aug 2016 16:35:16 -0700 Subject: [PATCH 4/5] Collapse redundancies --- source/simulator/initial_conditions.cc | 152 +++++++++---------------- 1 file changed, 52 insertions(+), 100 deletions(-) diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index f97f2319018..26ac8219354 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -53,102 +53,6 @@ namespace aspect initial_solution.reinit(system_rhs, false); -#if DEAL_II_VERSION_GTE(8,5,0) - - // Interpolate initial temperature - VectorFunctionFromScalarFunctionObject - temperature_function(std_cxx11::bind(&InitialConditions::Interface::initial_temperature, - std::ref(*initial_conditions), - std_cxx11::_1), - introspection.component_indices.temperature, - introspection.n_components); - - VectorTools::interpolate(*mapping, - dof_handler, - temperature_function, - initial_solution, - introspection.component_masks.temperature); - - // Interpolate initial compositions - for (unsigned int n=0; n - composition_function (std_cxx11::bind(&CompositionalInitialConditions::Interface::initial_composition, - std::ref(*compositional_initial_conditions), - std_cxx11::_1, - n), - introspection.component_indices.compositional_fields[n], - introspection.n_components); - - - VectorTools::interpolate(*mapping, - dof_handler, - composition_function, - initial_solution, - introspection.component_masks.compositional_fields[n]); - - const unsigned int base_element = advf.base_element(introspection); - - // get the composition support points - const std::vector > support_points - = finite_element.base_element(base_element).get_unit_support_points(); - Assert (support_points.size() != 0, - ExcInternalError()); - - // create an FEValues object to check compositions at support points - FEValues fe_values (*mapping, finite_element, - support_points, - update_quadrature_points); - - std::vector local_dof_indices (finite_element.dofs_per_cell); - - for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); - cell != dof_handler.end(); ++cell) - if (cell->is_locally_owned()) - { - fe_values.reinit (cell); - - // Go through the support points for each dof - for (unsigned int i=0; i0 && n == 0) - { - double sum = 0; - for (unsigned int m=0; minitial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]); - if (std::abs(sum) > 1.0+1e-6) - { - max_sum_comp = std::max(sum, max_sum_comp); - normalize_composition = true; - } - } - - } - } - - // if at least one processor decides that it needs - // to normalize, do the same on all processors. - if (Utilities::MPI::max (normalize_composition ? 1 : 0, - mpi_communicator) - == 1) - { - const double global_max - = Utilities::MPI::max (max_sum_comp, mpi_communicator); - - if (n==1) - pcout << "Sum of compositional fields is not one, fields will be normalized" - << std::endl; - - for (unsigned int m=0; m local_dof_indices (finite_element.dofs_per_cell); +#if DEAL_II_VERSION_GTE(8,5,0) + const VectorFunctionFromScalarFunctionObject &advf_init_function = + (advf.is_temperature()? + VectorFunctionFromScalarFunctionObject(std_cxx11::bind(&InitialConditions::Interface::initial_temperature, + std::ref(*initial_conditions), + std_cxx11::_1), + introspection.component_indices.temperature, + introspection.n_components): + VectorFunctionFromScalarFunctionObject(std_cxx11::bind(&CompositionalInitialConditions::Interface::initial_composition, + std::ref(*compositional_initial_conditions), + std_cxx11::_1, + n-1), + introspection.component_indices.compositional_fields[n-1], + introspection.n_components)); + + const ComponentMask advf_mask = + (advf.is_temperature()? + introspection.component_masks.temperature: + introspection.component_masks.compositional_fields[n-1]); + + VectorTools::interpolate(*mapping, + dof_handler, + advf_init_function, + initial_solution, + advf_mask); + + if (n == 1 && parameters.normalized_fields.size()>0) + for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + + // Go through the support points for each dof + for (unsigned int i=0; iinitial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]); + if (std::abs(sum) > 1.0+1e-6) + { + max_sum_comp = std::max(sum, max_sum_comp); + normalize_composition = true; + } + } + } +#else for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) if (cell->is_locally_owned()) @@ -226,6 +179,7 @@ namespace aspect } } +#endif initial_solution.compress(VectorOperation::insert); @@ -235,8 +189,6 @@ namespace aspect mpi_communicator) == 1) { - - const unsigned int blockidx = advf.block_index(introspection); const double global_max = Utilities::MPI::max (max_sum_comp, mpi_communicator); @@ -246,10 +198,10 @@ namespace aspect for (unsigned int m=0; m Date: Tue, 9 Aug 2016 10:17:52 -0700 Subject: [PATCH 5/5] Tidy up modification to initial_conditions.cc --- source/simulator/initial_conditions.cc | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index 26ac8219354..981d3341364 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -93,12 +93,14 @@ namespace aspect #if DEAL_II_VERSION_GTE(8,5,0) const VectorFunctionFromScalarFunctionObject &advf_init_function = - (advf.is_temperature()? + (advf.is_temperature() + ? VectorFunctionFromScalarFunctionObject(std_cxx11::bind(&InitialConditions::Interface::initial_temperature, std::ref(*initial_conditions), std_cxx11::_1), introspection.component_indices.temperature, - introspection.n_components): + introspection.n_components) + : VectorFunctionFromScalarFunctionObject(std_cxx11::bind(&CompositionalInitialConditions::Interface::initial_composition, std::ref(*compositional_initial_conditions), std_cxx11::_1, @@ -107,8 +109,10 @@ namespace aspect introspection.n_components)); const ComponentMask advf_mask = - (advf.is_temperature()? - introspection.component_masks.temperature: + (advf.is_temperature() + ? + introspection.component_masks.temperature + : introspection.component_masks.compositional_fields[n-1]); VectorTools::interpolate(*mapping, @@ -117,7 +121,7 @@ namespace aspect initial_solution, advf_mask); - if (n == 1 && parameters.normalized_fields.size()>0) + if (parameters.normalized_fields.size()>0 && n==1) for (typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell) if (cell->is_locally_owned()) @@ -131,8 +135,9 @@ namespace aspect // must not exceed one, this should be checked double sum = 0; for (unsigned int m=0; minitial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]); - if (std::abs(sum) > 1.0+1e-6) + sum += compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i), + parameters.normalized_fields[m]); + if (std::abs(sum) > 1.0+std::numeric_limits::epsilon()) { max_sum_comp = std::max(sum, max_sum_comp); normalize_composition = true; @@ -169,8 +174,9 @@ namespace aspect { double sum = 0; for (unsigned int m=0; minitial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]); - if (std::abs(sum) > 1.0+1e-6) + sum += compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i), + parameters.normalized_fields[m]); + if (std::abs(sum) > 1.0+std::numeric_limits::epsilon()) { max_sum_comp = std::max(sum, max_sum_comp); normalize_composition = true;