[cig-commits] commit 2407 by heister to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Sun Apr 6 15:02:04 PDT 2014
Revision 2407
ground work for nullspace removal: net rotation is working
U trunk/aspect/include/aspect/simulator.h
U trunk/aspect/source/simulator/core.cc
U trunk/aspect/source/simulator/helper_functions.cc
U trunk/aspect/source/simulator/parameters.cc
U trunk/aspect/source/simulator/solver.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2407&peg=2407
Diff:
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2014-04-04 03:40:38 UTC (rev 2406)
+++ trunk/aspect/include/aspect/simulator.h 2014-04-06 22:02:00 UTC (rev 2407)
@@ -38,6 +38,7 @@
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/mapping_q.h>
+#include <deal.II/base/tensor_function.h>
#include <aspect/global.h>
#include <aspect/simulator_access.h>
@@ -106,6 +107,18 @@
};
};
+ struct NullspaceRemoval
+ {
+ enum Kind
+ {
+ none = 0,
+ net_rotation = 0x1,
+ net_translation = 0x2,
+ angular_momentum = 0x4,
+ translational_momentum = 0x8
+ };
+ };
+
public:
/**
* A structure that holds run-time parameters. These parameters are all
@@ -197,6 +210,11 @@
*/
std::map<types::boundary_id, std::pair<std::string,std::string> > prescribed_velocity_boundary_indicators;
/**
+ * Selection of operations to perform to remove nullspace from
+ * velocity field.
+ */
+ typename NullspaceRemoval::Kind nullspace_removal;
+ /**
* @}
*/
@@ -899,7 +917,24 @@
*/
void denormalize_pressure(LinearAlgebra::BlockVector &vector);
+
/**
+ * interpolates the given function onto the velocity FE space and write it into the given vector.
+ */
+ void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
+ LinearAlgebra::Vector &vec);
+
+ /**
+ * Aets up data structures for null space removal. Called after every mesh refinement.
+ */
+ void setup_nullspace_removal();
+
+ /**
+ * Eliminate the nullspace of the velocity in the given vector @p vector.
+ */
+ void remove_nullspace(LinearAlgebra::BlockVector &vector);
+
+ /**
* Compute the maximal velocity throughout the domain. This is needed
* to compute the size of the time step.
*
@@ -1206,6 +1241,8 @@
bool rebuild_stokes_matrix;
bool rebuild_stokes_preconditioner;
+
+ std::vector<LinearAlgebra::Vector> net_rotations_translations;
/**
* @}
*/
Modified: trunk/aspect/source/simulator/core.cc
===================================================================
--- trunk/aspect/source/simulator/core.cc 2014-04-04 03:40:38 UTC (rev 2406)
+++ trunk/aspect/source/simulator/core.cc 2014-04-06 22:02:00 UTC (rev 2407)
@@ -1035,6 +1035,8 @@
old_solution = old_distributed_system;
}
+ setup_nullspace_removal();
+
computing_timer.exit_section();
}
Modified: trunk/aspect/source/simulator/helper_functions.cc
===================================================================
--- trunk/aspect/source/simulator/helper_functions.cc 2014-04-04 03:40:38 UTC (rev 2406)
+++ trunk/aspect/source/simulator/helper_functions.cc 2014-04-06 22:02:00 UTC (rev 2407)
@@ -469,7 +469,41 @@
}
+ template <int dim>
+ void Simulator<dim>::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
+ LinearAlgebra::Vector &vec)
+ {
+ ConstraintMatrix hanging_constraints(introspection.index_sets.system_relevant_set);
+ DoFTools::make_hanging_node_constraints(dof_handler, hanging_constraints);
+ hanging_constraints.close();
+ Assert(introspection.block_indices.velocities == 0, ExcNotImplemented());
+ const std::vector<Point<dim> > mesh_support_points = finite_element.base_element(0).get_unit_support_points();
+ FEValues<dim> mesh_points (mapping, finite_element, mesh_support_points, update_quadrature_points);
+ std::vector<unsigned int> cell_dof_indices (finite_element.dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for(; cell != endc; ++cell)
+ if(cell->is_locally_owned())
+ {
+ mesh_points.reinit(cell);
+ cell->get_dof_indices (cell_dof_indices);
+ for(unsigned int j=0; j<finite_element.base_element(0).dofs_per_cell; ++j)
+ for(unsigned int dir=0; dir<dim; ++dir)
+ {
+ unsigned int support_point_index
+ = finite_element.component_to_system_index(/*velocity component=*/ introspection.component_indices.velocities[dir],
+ /*dof index within component=*/ j);
+ Assert(introspection.block_indices.velocities == 0, ExcNotImplemented());
+ vec[cell_dof_indices[support_point_index]] = func.value(mesh_points.quadrature_point(j))[dir];
+ }
+ }
+
+ vec.compress(VectorOperation::insert);
+ hanging_constraints.distribute(vec);
+ }
+
/*
* normalize the pressure by calculating the surface integral of the pressure on the outer
* shell and subtracting this from all pressure nodes.
@@ -1038,7 +1072,8 @@
template void Simulator<dim>::compute_depth_average_Vp(std::vector<double> &values) const; \
template void Simulator<dim>::output_program_stats(); \
template void Simulator<dim>::output_statistics(); \
- template bool Simulator<dim>::stokes_matrix_depends_on_solution() const;
+ template bool Simulator<dim>::stokes_matrix_depends_on_solution() const; \
+ template void Simulator<dim>::interpolate_onto_velocity_system(const TensorFunction<1,dim> &func, LinearAlgebra::Vector &vec);
ASPECT_INSTANTIATE(INSTANTIATE)
}
Modified: trunk/aspect/source/simulator/parameters.cc
===================================================================
--- trunk/aspect/source/simulator/parameters.cc 2014-04-04 03:40:38 UTC (rev 2406)
+++ trunk/aspect/source/simulator/parameters.cc 2014-04-06 22:02:00 UTC (rev 2407)
@@ -342,6 +342,15 @@
"part of the boundary on which the velocity is to be zero with "
"the parameter ``Zero velocity boundary indicator'' in the "
"current parameter section.");
+
+ prm.declare_entry ("Remove nullspace", "",
+ Patterns::MultipleSelection("net rotation|net translation|angular momentum|translational momentum"),
+ "A selection of operations to remove certain parts of the nullspace from "
+ "the velocity after solving. For some geometries and certain boundary conditions "
+ "the velocity field is not uniquely determined but contains free translations "
+ "and or rotations.
"
+ "Note that while more than operation can be selected it only makes sense to "
+ "pick one rotational and one translational operation.");
}
prm.leave_subsection ();
@@ -700,6 +709,30 @@
prescribed_velocity_boundary_indicators[boundary_id] =
std::pair<std::string,std::string>(comp,value);
}
+
+ {
+ nullspace_removal = NullspaceRemoval::none;
+ std::vector<std::string> nullspace_names =
+ Utilities::split_string_list(prm.get("Remove nullspace"));
+ for (unsigned int i=0;i<nullspace_names.size();++i)
+ {
+ if (nullspace_names[i]=="net rotation")
+ nullspace_removal = typename NullspaceRemoval::Kind(
+ nullspace_removal | NullspaceRemoval::net_rotation);
+ else if (nullspace_names[i]=="net translation")
+ nullspace_removal = typename NullspaceRemoval::Kind(
+ nullspace_removal | NullspaceRemoval::net_translation);
+ else if (nullspace_names[i]=="angular momentum")
+ nullspace_removal = typename NullspaceRemoval::Kind(
+ nullspace_removal | NullspaceRemoval::angular_momentum);
+ else if (nullspace_names[i]=="translational momentum")
+ nullspace_removal = typename NullspaceRemoval::Kind(
+ nullspace_removal | NullspaceRemoval::translational_momentum);
+ else
+ AssertThrow(false, ExcInternalError());
+ }
+ }
+
}
prm.leave_subsection ();
Modified: trunk/aspect/source/simulator/solver.cc
===================================================================
--- trunk/aspect/source/simulator/solver.cc 2014-04-04 03:40:38 UTC (rev 2406)
+++ trunk/aspect/source/simulator/solver.cc 2014-04-06 22:02:00 UTC (rev 2407)
@@ -479,6 +479,8 @@
// now rescale the pressure back to real physical units
distributed_stokes_solution.block(1) *= pressure_scaling;
+ remove_nullspace(distributed_stokes_solution);
+
// then copy back the solution from the temporary (non-ghosted) vector
// into the ghosted one with all solution components
solution.block(0) = distributed_stokes_solution.block(0);
More information about the CIG-COMMITS
mailing list