diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h index 95a2cad44a3..7496ef68690 100644 --- a/include/aspect/simulator.h +++ b/include/aspect/simulator.h @@ -2097,7 +2097,6 @@ namespace aspect LinearAlgebra::BlockVector solution; LinearAlgebra::BlockVector old_solution; LinearAlgebra::BlockVector old_old_solution; - LinearAlgebra::BlockVector unlimited_solution; LinearAlgebra::BlockVector system_rhs; LinearAlgebra::BlockVector current_linearization_point; diff --git a/source/simulator/checkpoint_restart.cc b/source/simulator/checkpoint_restart.cc index 7f382d88f46..69a1e70281e 100644 --- a/source/simulator/checkpoint_restart.cc +++ b/source/simulator/checkpoint_restart.cc @@ -278,7 +278,7 @@ namespace aspect // save Triangulation and Solution vectors: { std::vector x_system - = { &solution, &old_solution, &old_old_solution, &unlimited_solution }; + = { &solution, &old_solution, &old_old_solution }; // If we are using a deforming mesh, include the mesh velocity, which uses the system dof handler if (parameters.mesh_deformation_enabled) @@ -540,11 +540,10 @@ namespace aspect LinearAlgebra::BlockVector distributed_system (system_rhs); LinearAlgebra::BlockVector old_distributed_system (system_rhs); LinearAlgebra::BlockVector old_old_distributed_system (system_rhs); - LinearAlgebra::BlockVector unlimited_distributed_system (system_rhs); LinearAlgebra::BlockVector distributed_mesh_velocity (system_rhs); std::vector x_system - = { &distributed_system, &old_distributed_system, &old_old_distributed_system, &unlimited_distributed_system }; + = { &distributed_system, &old_distributed_system, &old_old_distributed_system }; // If necessary, also include the mesh velocity for deserialization // with the system dof handler @@ -559,7 +558,6 @@ namespace aspect solution = distributed_system; old_solution = old_distributed_system; old_old_solution = old_old_distributed_system; - unlimited_solution = unlimited_distributed_system; if (parameters.mesh_deformation_enabled) { diff --git a/source/simulator/core.cc b/source/simulator/core.cc index eab5f37ff0e..e2ce9bfd73d 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -1486,7 +1486,6 @@ namespace aspect solution.reinit(introspection.index_sets.system_partitioning, introspection.index_sets.system_relevant_partitioning, mpi_communicator); old_solution.reinit(introspection.index_sets.system_partitioning, introspection.index_sets.system_relevant_partitioning, mpi_communicator); old_old_solution.reinit(introspection.index_sets.system_partitioning, introspection.index_sets.system_relevant_partitioning, mpi_communicator); - unlimited_solution.reinit(introspection.index_sets.system_partitioning, introspection.index_sets.system_relevant_partitioning, mpi_communicator); current_linearization_point.reinit (introspection.index_sets.system_partitioning, introspection.index_sets.system_relevant_partitioning, mpi_communicator); if (parameters.use_operator_splitting) @@ -1670,7 +1669,7 @@ namespace aspect // Next set up whatever is necessary to transfer the solution from old // to new mesh: std::vector x_system - = { &solution, &old_solution, &unlimited_solution }; + = { &solution, &old_solution }; if (parameters.mesh_deformation_enabled) x_system.push_back(&mesh_deformation->mesh_velocity); @@ -1728,17 +1727,15 @@ namespace aspect LinearAlgebra::BlockVector distributed_system; LinearAlgebra::BlockVector old_distributed_system; - LinearAlgebra::BlockVector unlimited_distributed_system; LinearAlgebra::BlockVector distributed_mesh_velocity; distributed_system.reinit(introspection.index_sets.system_partitioning, mpi_communicator); old_distributed_system.reinit(introspection.index_sets.system_partitioning, mpi_communicator); - unlimited_distributed_system.reinit(introspection.index_sets.system_partitioning, mpi_communicator); if (parameters.mesh_deformation_enabled) distributed_mesh_velocity.reinit(introspection.index_sets.system_partitioning, mpi_communicator); std::vector system_tmp - = { &distributed_system, &old_distributed_system, &unlimited_distributed_system}; + = { &distributed_system, &old_distributed_system}; if (parameters.mesh_deformation_enabled) system_tmp.push_back(&distributed_mesh_velocity); @@ -1765,9 +1762,6 @@ namespace aspect constraints.distribute (old_distributed_system); old_solution = old_distributed_system; - constraints.distribute (unlimited_distributed_system); - unlimited_solution = unlimited_distributed_system; - // We need the current linearization point at the start of the new time step // when we set the boundary conditions for advected fields (to determine parts // of the boundary with outflow). Therefore, we here set it to the solution diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index ba48d434534..91bc59bb71b 100644 --- a/source/simulator/helper_functions.cc +++ b/source/simulator/helper_functions.cc @@ -1618,17 +1618,6 @@ namespace aspect { solution.block(block_index) = distributed_vector.block(block_index); - // If a limiter is used for a Discontinuous Galerkin field, - // also update the stored unlimited solution. - // TODO the update is computed using the limited solution, - // is that consistent? - if ((parameters.use_discontinuous_temperature_discretization - && parameters.use_limiter_for_discontinuous_temperature_solution) - || - (parameters.use_discontinuous_composition_discretization - && parameters.use_limiter_for_discontinuous_composition_solution)) - unlimited_solution.block(block_index) = distributed_vector.block(block_index); - // we have to update the old solution with our reaction update too // so that the advection scheme will have the correct time stepping in the next step LinearAlgebra::BlockVector tmp; diff --git a/source/simulator/initial_conditions.cc b/source/simulator/initial_conditions.cc index a2d60a8a384..e1733dc2192 100644 --- a/source/simulator/initial_conditions.cc +++ b/source/simulator/initial_conditions.cc @@ -242,7 +242,6 @@ namespace aspect const unsigned int blockidx = advf.block_index(introspection); solution.block(blockidx) = initial_solution.block(blockidx); - unlimited_solution.block(blockidx) = initial_solution.block(blockidx); old_solution.block(blockidx) = initial_solution.block(blockidx); old_old_solution.block(blockidx) = initial_solution.block(blockidx); current_linearization_point.block(blockidx) = initial_solution.block(blockidx); diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc index e6702d68b2d..ed35b91674d 100644 --- a/source/simulator/solver.cc +++ b/source/simulator/solver.cc @@ -502,47 +502,10 @@ namespace aspect // Compute the residual before we solve and return this at the end. // This is used in the nonlinear solver. - double initial_residual = system_matrix.block(block_idx,block_idx).residual - (temp, - distributed_solution.block(block_idx), - system_rhs.block(block_idx)); - - // If a limiter is used for a Discontinuous Galerkin field, - // use the unlimited solution to compute the residual. - // Only do this after the first nonlinear iteration, - // because unlike the current_linearization_point, - // the unlimited solution is not extrapolated from the two - // previous timesteps. After solving once within a timestep, - // we have an unlimited solution we can use. - if (((advection_field.is_temperature() - && parameters.use_discontinuous_temperature_discretization - && parameters.use_limiter_for_discontinuous_temperature_solution) - || - (!advection_field.is_temperature() - && parameters.use_discontinuous_composition_discretization - && parameters.use_limiter_for_discontinuous_composition_solution)) - && nonlinear_iteration > 0) - { - LinearAlgebra::BlockVector distributed_unlimited_solution ( - introspection.index_sets.system_partitioning, - mpi_communicator); - - distributed_unlimited_solution.block(block_idx) = unlimited_solution.block (block_idx); - - // Temporary vector to hold the residual, we don't need a BlockVector here. - LinearAlgebra::Vector temp2 ( - introspection.index_sets.system_partitioning[block_idx], - mpi_communicator); - - current_constraints.set_zero(distributed_unlimited_solution); - - const double unlimited_initial_residual = system_matrix.block(block_idx,block_idx).residual - (temp2, - distributed_unlimited_solution.block(block_idx), - system_rhs.block(block_idx)); - - initial_residual = unlimited_initial_residual; - } + const double initial_residual = system_matrix.block(block_idx,block_idx).residual + (temp, + distributed_solution.block(block_idx), + system_rhs.block(block_idx)); // solve the linear system: try