Dear Wolfgang,
I attach a minimal working program (sorry its still quite lengthy - a lot of it
is outputting the solution and includes etc.). I have only run this in 2D at
the moment (its currently set to run in 2D), so I'm not sure that it is working
properly in 3D - just in case!
The BlockVector theta_solution consists of theta_n and theta_s, and it is
theta_s that is causing the problem. This should be either 0.03 or 0 (set at
the top by ThetaInitialValues in EquationData and accessed through
VectorTools::project in the "run" part of the program). At the moment
BlockVector stokes_solution does nothing, but will contain 2 velocity fields
and pressure, I have left it intact as I don't know whether it has something to
do with the overshoots or not.
Where the over-shoots occur also depends on the initial refinement (there is
only initial refinement in this program, which occurs before the values are
projected onto the grid.)
Any help would be greatly appreciated.
Many Thanks,
Katie Leonard
DPhil student in Computational Biology,
The University of Oxford.
________________________________________
From: Wolfgang Bangerth [[email protected]]
Sent: 09 November 2011 16:08
To: Katie Leonard
Cc: Guido Kanschat; [email protected]
Subject: Re: [deal.II] Projecting Initial Values without interpolation
> Sorry for such a ridiculously long delay. I attach a picture to
> demonstrate what I mean.
>
> This picture should show $\theta_s$ either taking a value of 0.03 or 0,
> depending on spatial co-ordinates. I don't know how values greater than
> 0.03 are appearing, as I'm doing nothing but projecting the values onto an
> already refined grid.
I would say write the minimal program that shows this problem (should be
doable in one or two pages of code) and send it to us. It's indeed strange
that there are only over- but not undershoots and that the overshoots only
happen on one side of the square...
W.
-------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.tamu.edu/~bangerth/
// K H L Leonard modified from step-31.cc
#include <base/quadrature_lib.h>
#include <base/logstream.h>
#include <base/utilities.h>
#include <lac/full_matrix.h>
#include <lac/solver_gmres.h>
#include <lac/solver_cg.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/constraint_matrix.h>
#include <lac/precondition.h>
#include <grid/tria.h>
#include <grid/grid_generator.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <grid/tria_boundary_lib.h>
#include <grid/grid_tools.h>
#include <grid/grid_refinement.h>
#include <dofs/dof_handler.h>
#include <dofs/dof_renumbering.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
#include <fe/fe_q.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <numerics/vectors.h>
#include <numerics/data_out.h>
#include <numerics/error_estimator.h>
#include <numerics/solution_transfer.h>
#include <base/parameter_handler.h>
#include <base/tensor_function.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_sparse_matrix.h>
#include <lac/trilinos_vector.h>
#include <lac/trilinos_block_vector.h>
#include <lac/trilinos_precondition.h>
#include <fstream>
#include <sstream>
#include <iostream>
using namespace dealii;
namespace EquationData
{
const double IC_scale = 50; // parameter controlling steepness of
switch-like initial conditions
double radius_scaffold = 0.2;
double radius_cell_IC = 0.2;
double height_culture = 0.5; // should be given in parameter file....2*"x
axis limits"
double height_scaffold = 0.2;
template <int dim>
class ThetaInitialValues : public Function<dim>
{
public:
ThetaInitialValues () : Function<dim>(2) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
ThetaInitialValues<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
if(component == 0) // theta_n
{
return 0.1*(tanh(IC_scale*(p[0]+0.5*height_culture)) - tanh(IC_scale*(p[0]
- height_scaffold + 0.5*height_culture)) )*(tanh(IC_scale*(p[1]
+radius_cell_IC)) - tanh(IC_scale*(p[1] - radius_cell_IC)));// + 0.000001;
}
else if(component == 1) // theta_s
{
if(p[1]>=-radius_scaffold && p[1] <=radius_scaffold && p[0] <=
height_scaffold-0.5*height_culture)
{
return 0.03;
}
else
{
return 0;
}
}
else
{
return 0;
}
}
template <int dim>
void
ThetaInitialValues<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = ThetaInitialValues<dim>::value (p, c);
}
}
template <int dim>
class BoussinesqFlowProblem
{
public:
BoussinesqFlowProblem ();
void run ();
private:
void setup_dofs ();
void output_results (const unsigned int initial_refinement) const;
Triangulation<dim> triangulation;
const unsigned int stokes_degree;
FESystem<dim> stokes_fe;
DoFHandler<dim> stokes_dof_handler;
ConstraintMatrix stokes_constraints;
std::vector<unsigned int> stokes_block_sizes;
TrilinosWrappers::BlockVector stokes_solution;
const unsigned int theta_degree;
FESystem<dim> theta_fe;
DoFHandler<dim> theta_dof_handler;
ConstraintMatrix theta_constraints;
std::vector<unsigned int> theta_block_sizes;
TrilinosWrappers::BlockVector theta_solution;
};
template <int dim>
BoussinesqFlowProblem<dim>::BoussinesqFlowProblem ()
:
triangulation (Triangulation<dim>::maximum_smoothing),
stokes_degree (1),
stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
FE_Q<dim>(stokes_degree+1), dim,
FE_Q<dim>(stokes_degree), 1),
stokes_dof_handler (triangulation),
theta_degree (2),
theta_fe (FE_Q<dim>(theta_degree), 1,
FE_Q<dim>(theta_degree), 1),
theta_dof_handler (triangulation)
{}
template <int dim>
void BoussinesqFlowProblem<dim>::setup_dofs ()
{
std::vector<unsigned int> stokes_sub_blocks (2*dim+1,0); // 2 velocity
systems and pressure
stokes_sub_blocks[2*dim] = 1; //pressure component
{
stokes_dof_handler.distribute_dofs (stokes_fe);
DoFRenumbering::component_wise (stokes_dof_handler, stokes_sub_blocks);
{
stokes_constraints.clear ();
DoFTools::make_hanging_node_constraints (stokes_dof_handler,
stokes_constraints);
}
stokes_constraints.close ();
}
std::vector<unsigned int> theta_sub_blocks (2,0);
theta_sub_blocks[1] = 1;
{
theta_dof_handler.distribute_dofs (theta_fe);
DoFRenumbering::component_wise (theta_dof_handler, theta_sub_blocks);
theta_constraints.clear ();
DoFTools::make_hanging_node_constraints (theta_dof_handler,
theta_constraints);
theta_constraints.close ();
}
std::vector<unsigned int> stokes_dofs_per_block (2);
DoFTools::count_dofs_per_block (stokes_dof_handler, stokes_dofs_per_block,
stokes_sub_blocks);
std::vector<unsigned int> theta_dofs_per_block (2);
DoFTools::count_dofs_per_block (theta_dof_handler, theta_dofs_per_block,
theta_sub_blocks);
const unsigned int n_u = stokes_dofs_per_block[0],
n_p = stokes_dofs_per_block[1],
n_theta_n = theta_dofs_per_block[0],
n_theta_s = theta_dofs_per_block[1];
stokes_block_sizes.resize (2);
stokes_block_sizes[0] = n_u;
stokes_block_sizes[1] = n_p;
theta_block_sizes.resize (2);
theta_block_sizes[0] = n_theta_n;
theta_block_sizes[1] = n_theta_s;
stokes_solution.reinit (stokes_block_sizes);
theta_solution.reinit (theta_block_sizes);
}
template <int dim>
void BoussinesqFlowProblem<dim>::output_results (const unsigned int
initial_refinement) const
{
const FESystem<dim> joint_fe (stokes_fe, 1,
theta_fe, 1);
DoFHandler<dim> joint_dof_handler (triangulation);
joint_dof_handler.distribute_dofs (joint_fe);
Assert (joint_dof_handler.n_dofs() ==
stokes_dof_handler.n_dofs() + theta_dof_handler.n_dofs(),
ExcInternalError());
Vector<double> joint_solution (joint_dof_handler.n_dofs());
{
std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
std::vector<unsigned int> local_stokes_dof_indices
(stokes_fe.dofs_per_cell);
std::vector<unsigned int> local_theta_dof_indices (theta_fe.dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
joint_cell = joint_dof_handler.begin_active(),
joint_endc = joint_dof_handler.end(),
stokes_cell = stokes_dof_handler.begin_active(),
theta_cell = theta_dof_handler.begin_active();
for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++theta_cell)
{
joint_cell->get_dof_indices (local_joint_dof_indices);
stokes_cell->get_dof_indices (local_stokes_dof_indices);
theta_cell->get_dof_indices (local_theta_dof_indices);
for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
if (joint_fe.system_to_base_index(i).first.first == 0)
{
Assert (joint_fe.system_to_base_index(i).second
<
local_stokes_dof_indices.size(),
ExcInternalError());
joint_solution(local_joint_dof_indices[i])
=
stokes_solution(local_stokes_dof_indices[joint_fe.system_to_base_index(i).second]);
}
else
{
Assert (joint_fe.system_to_base_index(i).first.first == 1,
ExcInternalError());
Assert (joint_fe.system_to_base_index(i).second
<
local_theta_dof_indices.size(),
ExcInternalError());
joint_solution(local_joint_dof_indices[i])
=
theta_solution(local_theta_dof_indices[joint_fe.system_to_base_index(i).second]);
}
}
}
std::vector<std::string> joint_solution_names;
switch (dim)
{
case 2:
joint_solution_names.push_back("u_n");
joint_solution_names.push_back("v_n");
joint_solution_names.push_back("u_w");
joint_solution_names.push_back("v_w");
joint_solution_names.push_back("p");
joint_solution_names.push_back("theta_n");
joint_solution_names.push_back("theta_s");
break;
case 3:
joint_solution_names.push_back("u_n");
joint_solution_names.push_back("v_n");
joint_solution_names.push_back("w_n");
joint_solution_names.push_back("u_w");
joint_solution_names.push_back("v_w");
joint_solution_names.push_back("w_w");
joint_solution_names.push_back("p");
joint_solution_names.push_back("theta_n");
joint_solution_names.push_back("theta_s");
break;
default:
Assert (false, ExcNotImplemented());
}
DataOut<dim> data_out;
data_out.attach_dof_handler (joint_dof_handler);
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation
(2*dim+3, DataComponentInterpretation::component_is_scalar);
for (unsigned int i=0; i<2*dim; ++i)
data_component_interpretation[i]
= DataComponentInterpretation::component_is_part_of_vector;
data_out.add_data_vector (joint_solution, joint_solution_names,
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.build_patches (std::min(stokes_degree, theta_degree));
std::ostringstream filename;
if(dim == 2)
{
filename << "modified_CPDG_hydrostatic_simple_"
<< Utilities::int_to_string(EquationData::IC_scale,2)
<< "_"
<< Utilities::int_to_string(initial_refinement, 1)
<< ".gpl";
std::ofstream output (filename.str().c_str());
data_out.write_gnuplot(output);
}
else
{
filename << "hydrostatic_simple-"
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
}
template <int dim>
void BoussinesqFlowProblem<dim>::run ()
{
const unsigned int initial_refinement = (dim == 2 ? 7: 2);
const double radius = 0.5;
const double xvalue = 2;
GridGenerator::cylinder(triangulation, radius, xvalue);
for (typename Triangulation<dim>::active_cell_iterator
cell=triangulation.begin_active();
cell!=triangulation.end(); ++cell)
{
if(dim == 3)
{
static const CylinderBoundary<dim> cylinder_boundary (radius, 0);
triangulation.set_boundary (0, cylinder_boundary);
for (typename Triangulation<dim>::active_face_iterator
face = triangulation.begin_active_face();
face!=triangulation.end_face(); ++face)
if(face->at_boundary())
if(face->boundary_indicator() == 0)
for(unsigned int edge = 0;
edge<GeometryInfo<dim>::lines_per_face; ++edge)
face->line(edge)->set_boundary_indicator(face->boundary_indicator());
}
}
triangulation.refine_global (initial_refinement);
setup_dofs();
unsigned int pre_refinement_step = 0;
VectorTools::project (theta_dof_handler,
theta_constraints,
QGauss<dim>(theta_degree+2),
EquationData::ThetaInitialValues<dim>(),
theta_solution);
output_results (initial_refinement);
}
int main (int argc, char *argv[])
{
try
{
deallog.depth_console (0);
Utilities::System::MPI_InitFinalize mpi_initialization (argc, argv);
BoussinesqFlowProblem<2> flow_problem;
flow_problem.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
