Hi folks,
I just started working with deal.II and do experience some problems with
solving the following equations:
\nabla n = 0
\nabla \vec{E} = -e(n - \Theta(x))
For simplicity, I started off in 1d where I can solve the equations
analytically. The problem I'm stuck with is
"Exception on processing:
Iterative method reported convergence failure in step 1 with residual nan"
I attached the program so you can have a look. What I actually did is to
take the Step-by-Step example on periodic boundary conditions (since in
higher dimensions I want to apply them) and modified it using a FESystem
and of course inserted my equations and left out the periodics as I
don't need them in 1d.
I tried using different Solver classes ending up with the same error.
So, I guess there's something wrong with setting up the system, but
after rereading the examples and the explanations on the module for
solving vector-valued eqns, I have no idea what's going wrong.
Am I missing something?
I would appreciate any help on this. Thanks.
Timo
//============================================================================
// Name : HDSSolver1d.cpp
// Author : koellner
// Version :
// Copyright :
// Description : Hello World in C++, Ansi-style
//============================================================================
#include "dealii/grid/tria.h"
#include "dealii/grid/grid_generator.h"
#include "dealii/grid/grid_out.h"
#include "dealii/grid/tria_boundary.h"
#include "dealii/fe/fe_q.h"
#include "dealii/fe/fe_system.h"
#include "dealii/fe/fe_values.h"
#include "dealii/dofs/dof_handler.h"
#include "dealii/lac/sparse_matrix.h"
#include "dealii/lac/sparsity_pattern.h"
#include "dealii/lac/constraint_matrix.h"
#include "dealii/lac/vector.h"
#include "dealii/lac/compressed_sparsity_pattern.h"
#include "dealii/lac/full_matrix.h"
#include "dealii/lac/precondition.h"
#include "dealii/lac/solver_cg.h"
#include "dealii/lac/solver_gmres.h"
#include "dealii/lac/solver_minres.h"
#include "dealii/lac/solver_control.h"
#include "dealii/numerics/matrices.h"
#include "dealii/numerics/vectors.h"
#include "dealii/dofs/dof_tools.h"
#include "dealii/base/quadrature.h"
#include "dealii/base/quadrature_lib.h"
#include "dealii/numerics/data_out.h"
#include "ThetaFunction.h"
#include <iostream>
#include <fstream>
#include <math.h>
class HDSProblem
{
public:
HDSProblem();
void run();
private:
void makeGridAndDofs();
void assembleSystem();
void solve();
void outputResults();
dealii::Triangulation<1> tria;
dealii::FESystem<1> fe;
dealii::DoFHandler<1> dofHandler;
dealii::SparsityPattern sPattern;
dealii::SparseMatrix<double> systemMatrix;
dealii::ConstraintMatrix constraints;
dealii::Vector<double> solution;
dealii::Vector<double> rhs;
double e;
double nbar;
double hbar;
double mstar;
double pi;
double xi;
double alpha;
double beta;
};
HDSProblem::HDSProblem()
: fe(dealii::FE_Q<1>(2), 2),
dofHandler(tria)
{
e = 0.0180951262047;
nbar = 58500000000;
hbar = 0.00357065584419368872;
mstar = 8875.05148793008989951356;
pi = 3.14159265358979323846;
xi = pow(3*pi*pi,2/3)/(5*mstar)*hbar*hbar;
alpha = (-3*e)/(pow(5*xi*nbar,2/3));
beta = -e*nbar;
}
void HDSProblem::makeGridAndDofs()
{
std::cout << "Making grid and dofs..." << std::endl;
dealii::GridGenerator::hyper_cube(tria, -25e-4, 25e-4);
//set surface flags to determine which to apply dirichlet and periodic
tria.begin_active()->face(0)->set_boundary_indicator(1); //left surface
tria.begin_active()->face(1)->set_boundary_indicator(1); //right surface
tria.begin_active()->face(2)->set_boundary_indicator(2); //bottom surface
tria.begin_active()->face(3)->set_boundary_indicator(2); //top
//refine grid
tria.refine_global(3);
dofHandler.distribute_dofs(fe);
std::cout << "Number of active cells: "
<< tria.n_active_cells ()
<< std::endl
<< "Dofs: " << dofHandler.n_dofs ()
<< std::endl;
constraints.clear();
//dirichlet boundaries on left and right surface
std::vector<bool> comps(2,true);
comps[1] = false;
dealii::VectorTools::interpolate_boundary_values(dofHandler, 1,
dealii::ZeroFunction<1>(2),
constraints,
comps);
constraints.close();
dealii::CompressedSparsityPattern cSparsity(dofHandler.n_dofs(), dofHandler.n_dofs());
dealii::DoFTools::make_sparsity_pattern(dofHandler, cSparsity, constraints, false);
cSparsity.compress();
sPattern.copy_from(cSparsity);
systemMatrix.reinit(sPattern);
solution.reinit(dofHandler.n_dofs());
rhs.reinit(dofHandler.n_dofs());
}
void HDSProblem::assembleSystem()
{
std::cout << "Assemble System..." << std::endl;
const dealii::FEValuesExtractors::Vector electricField(0);
const dealii::FEValuesExtractors::Scalar density(1);
dealii::QGauss<1> quad(2);
dealii::FEValues<1> feValues(fe, quad, dealii::update_JxW_values | dealii::update_gradients |
dealii::update_values | dealii::update_quadrature_points);
const unsigned int nDofsPerCell = fe.dofs_per_cell;
const unsigned int nQPoints = quad.size();
dealii::FullMatrix<double> cellMatrix(nDofsPerCell, nDofsPerCell);
dealii::Vector<double> cellRhs(nDofsPerCell);
std::vector<unsigned int> dofsLocalToGlobal(nDofsPerCell);
const ThetaFunction<1> tf;
dealii::DoFHandler<1>::active_cell_iterator cell = dofHandler.begin_active(),
end = dofHandler.end();
for(; cell != end; ++cell)
{
feValues.reinit(cell);
cellMatrix = 0;
cellRhs = 0;
for(unsigned int q = 0; q < nQPoints; ++q)
{
for(unsigned int i = 0; i < nDofsPerCell; ++i)
{
for(unsigned int j = 0; j < nDofsPerCell; ++j)
{
cellMatrix(i,j) += (feValues[electricField].value(i,q) * feValues[density].gradient(j,q) +
feValues[density].value(i,q) * feValues[electricField].divergence(j,q) +
feValues[density].value(i,q) * feValues[density].value(j,q) * beta) * feValues.JxW(q);
}
cellRhs(i) += (feValues[density].value(i,q) * beta * tf.value(feValues.quadrature_point(q))) * feValues.JxW(q);
}
cell->get_dof_indices(dofsLocalToGlobal);
constraints.distribute_local_to_global(cellMatrix, cellRhs, dofsLocalToGlobal, systemMatrix, rhs);
}
}
}
void HDSProblem::solve()
{
std::cout << "Solving system..." << std::endl;
dealii::SolverControl solverControl (dofHandler.n_dofs (), 1e-12);
solverControl.enable_history_data();
solverControl.log_history(true);
solverControl.log_frequency(1);
solverControl.log_result(true);
dealii::PreconditionSSOR<dealii::SparseMatrix<double> > precondition;
precondition.initialize (systemMatrix);
dealii::SolverCG<> cg (solverControl);
cg.solve (systemMatrix, solution, rhs, precondition);
constraints.distribute (solution);
std::cout << " " << solverControl.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;}
void HDSProblem::run()
{
makeGridAndDofs();
assembleSystem();
solve();
}
int main ()
{
try
{
dealii::deallog.depth_console (0);
HDSProblem hds;
hds.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
