Undo using unlimited solution for DG residual

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;

source/simulator/checkpoint_restart.cc

@@ -278,7 +278,7 @@ namespace aspect
     // save Triangulation and Solution vectors:
     {
       std::vector<const LinearAlgebra::BlockVector *> 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<LinearAlgebra::BlockVector *> 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)
       {

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<const LinearAlgebra::BlockVector *> 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<LinearAlgebra::BlockVector *> 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

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;

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);

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