Re: [deal.II] Re: Usage of the laplace-matrix in example 23
Dulcimer, it is difficult for any of us to just look at a piece of code and tell where exactly things are going wrong. You'll have to learn to debug what is happening in your case, and for this it is easiest if you make the problem as simple as possible -- for example, use zero boundary values, zero right hand sides, no hanging nodes, etc. You can then reduce the amount of code that deal with each of these cases. This reduces the number of possible sources of the problem. From a cursory look at your code, I see that you are adding contributions to the rhs vectors and the system matrix on each cell. But I don't see that you set these vectors and the matrix to zero before your loop over all cells. I may be missing this, but it seems like a bug. Best W. -- Wolfgang Bangerth email: bange...@colostate.edu www: http://www.math.colostate.edu/~bangerth/ -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] Re: Usage of the laplace-matrix in example 23
Please excuse me for the previous post.
Here is the corrected function:
a) Added comments for what I have added and removed
b) Also while assembling the matrix_v. I was using the matrix_u (mistake
which arose from my copy paste technique.)
template
void WaveEquation::run ()
{
setup_system();
/*New*/
QGauss quadrature_formula(3);
/*New*/
VectorTools::project (dof_handler, constraints, QGauss(3),
InitialValuesU(),
old_solution_u);
VectorTools::project (dof_handler, constraints, QGauss(3),
InitialValuesV(),
old_solution_v);
Vector tmp (solution_u.size());
Vector forcing_terms (solution_u.size());
//***
// Code Substituted
//***
FEValues fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
//***
std::cout << "Dofs per cell : " << dofs_per_cell << std::endl;
std::cout << "Quadrature formula size : " << n_q_points << std::endl;
for (; time<=5; time+=time_step, ++timestep_number)
{
std::cout << "Time step " << timestep_number
<< " at t=" << time
<< std::endl;
//***
// Code Substituted
//***
FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
Vector cell_rhs(dofs_per_cell);
Vector cell_rhs_1(dofs_per_cell);
Vector cell_rhs_2(dofs_per_cell);
std::vector
local_dof_indices(dofs_per_cell);
std::vector old_u_q(n_q_points);
std::vector old_v_q(n_q_points);
for(const auto &cell: dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
cell_rhs_1 = 0;
cell_rhs_2 = 0;
fe_values.get_function_values(old_solution_u, old_u_q);
fe_values.get_function_values(old_solution_v, old_v_q);
for(unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
double old_u = old_u_q[q_index];
double old_v = old_v_q[q_index];
for(unsigned int i = 0; i < dofs_per_cell; ++i)
{
for(unsigned int j = 0; j < dofs_per_cell; ++j)
{
/*(M + k^2*theta^2*A)*/
cell_matrix(i, j) += (fe_values.shape_value(i,
q_index) *
fe_values.shape_value(j,
q_index)
+
time_step*time_step*
theta * theta *
fe_values.shape_grad(i,
q_index) *
fe_values.shape_grad(j,
q_index)) *
fe_values.JxW(q_index);
cell_rhs_1(i) +=
/*(M - k^2*theta*(1-theta)*A)U^(n-1)*/
(fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index)
-
time_step*time_step*
theta * (1 - theta)*
fe_values.shape_grad(i, q_index) *
fe_values.shape_grad(j, q_index)
) *
old_u;
cell_rhs_2(i) +=
/*k*M*V^(n-1)*/
time_step *
fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index) *
old_v;
}
cell_rhs(i) += (cell_rhs_1(i) + cell_rhs_2(i) *
fe_values.JxW (q_index));
}
}
cell->get_dof_indices (local_dof_indices);
for(unsigned int i = 0; i < dofs_per_cell; i++)
{
for (unsigned int j = 0; j < dofs_per_cell; j++)
{
matrix_u.add (lo
Re: [deal.II] Re: Usage of the laplace-matrix in example 23
Hello Professor,
I deleted my question since I thought I had found the answer(in a previous
post of yours) But I guess I am still stuck in the same problem.
What I have done in step 23 is to only make changes in the existing run
function and assemble my lhs matrices and rhs vectors cell wise.
Here is the run function I have tried to implement but the total energy
keeps increasing. Also, since the F values are considered zero, I thought I
should bother myself with F before I get this running. For the past three
days I have been trying to find what I have been missing but, alas :
Thanking you for your time.
template
void WaveEquation::run ()
{
setup_system();
QGauss quadrature_formula(3);
VectorTools::project (dof_handler, constraints, QGauss(3),
InitialValuesU(),
old_solution_u);
VectorTools::project (dof_handler, constraints, QGauss(3),
InitialValuesV(),
old_solution_v);
Vector tmp (solution_u.size());
Vector forcing_terms (solution_u.size());
//***
// Code Substituted
//***
FEValues fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
//***
std::cout << "Dofs per cell : " << dofs_per_cell << std::endl;
std::cout << "Quadrature formula size : " << n_q_points << std::endl;
for (; time<=5; time+=time_step, ++timestep_number)
{
std::cout << "Time step " << timestep_number
<< " at t=" << time
<< std::endl;
//***
// Code Substituted
//***
FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
Vector cell_rhs(dofs_per_cell);
Vector cell_rhs_1(dofs_per_cell);
Vector cell_rhs_2(dofs_per_cell);
std::vector
local_dof_indices(dofs_per_cell);
std::vector old_u_q(n_q_points);
std::vector old_v_q(n_q_points);
for(const auto &cell: dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
cell_rhs_1 = 0;
cell_rhs_2 = 0;
fe_values.get_function_values(old_solution_u, old_u_q);
fe_values.get_function_values(old_solution_v, old_v_q);
for(unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
double old_u = old_u_q[q_index];
double old_v = old_v_q[q_index];
for(unsigned int i = 0; i < dofs_per_cell; ++i)
{
for(unsigned int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix(i, j) += (fe_values.shape_value(i,
q_index) *
fe_values.shape_value(j,
q_index)
+
time_step*time_step*
theta * theta *
fe_values.shape_grad(i,
q_index) *
fe_values.shape_grad(j,
q_index)) *
fe_values.JxW(q_index);
cell_rhs_1(i) +=
/*(M - k^2*theta*(1-theta)*A)U^(n-1)*/
(fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index)
-
time_step*time_step*
theta * (1 - theta)*
fe_values.shape_grad(i, q_index) *
fe_values.shape_grad(j, q_index)
) *
old_u;
cell_rhs_2(i) +=
/*k*M*V^(n-1)*/
time_step *
fe_values.shape_value(i, q_index) *
fe_values.shape_value(j, q_index) *
old_v;
}
cell_rhs(i) += (cell_rhs_1(i) + cell_rhs_2(i) *
Re: [deal.II] Re: Usage of the laplace-matrix in example 23
On 01/23/2018 10:35 AM, Dulcimer0909 wrote: If I do go ahead and replace code, so that it does a cell by cell assembly, I am a bit lost on how I would store the old_solution (U^(n-1)) for each cell and retrieve it during the assembly for the Rhs. Dulcimer -- can you elaborate? It's not clear to me how the two are related. In the thread you comment on, the question is about how to assemble the matrices, but these matrices do not actually depend on the previous solution. In any case, the program of course stores the old solution and so any code you write will have access to it. So it's not entirely clear to me what you mean in your question :-( Best W. -- Wolfgang Bangerth email: bange...@colostate.edu www: http://www.math.colostate.edu/~bangerth/ -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Usage of the laplace-matrix in example 23
Hello all, An additional question regarding this thread: If I do go ahead and replace code, so that it does a cell by cell assembly, I am a bit lost on how I would store the old_solution (U^(n-1)) for each cell and retrieve it during the assembly for the Rhs. grateful if anyone could help. Dulcimer On Tuesday, August 22, 2017 at 4:50:03 PM UTC+2, Daniel Arndt wrote: > > Maxi, > > citing from the documentation of step-23[1]: > "In a very similar vein, we are also too lazy to write the code to > assemble mass and Laplace matrices, although it would have only taken > copying the relevant code from any number of previous tutorial programs. > Rather, we want to focus on the things that are truly new to this program > and therefore use the MatrixCreator::create_mass_matrix and > MatrixCreator::create_laplace_matrix functions." > It comes just handy that we have an easy way to get the required matrices > without having to write the assembly loops ourselves. > Mutiplying the finite element vector old_solution_u by the laplace matrix > is the same as assembling the right-hand side (\nabla phi, \nabla > old_solution_u). > > The point in this example is that you don't have to assemble the laplace > matrix for each time step anew. How you do this is totally up to you. In > create MatrixCreator::create_laplace_matrix we apply some more optimizations > so this might be faster then writing the assembly loop yourself. In the > end, you should have a lot of time steps and assembling this matrix once > should negligible in comparison to the time you spend in the solver. > Hence, my advice would be to use whatever is easier for you. > > Best, > Daniel > > [1] https://www.dealii.org/8.5.0/doxygen/deal.II/step_23.html#Includefiles > > > Am Freitag, 18. August 2017 17:26:37 UTC+2 schrieb Maxi Miller: >> >> I am trying to develop my own code based on example 23 (\partial_t U = >> i/(2k)\nabla^2U), and am trying to understand the purpose of the >> laplace-matrix. Why is it used in order to multiply in both sides, in >> comparison with the weak formulation (\nabla\phi_i,\nabla\phi_j) and the >> assembly of the usual cell matrix? Assumed I have to split my U into a real >> and an imaginary value, is the approach with the Laplace-matrix better, or >> the split into the weak formulation and the accompanying matrices on both >> sides? >> > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Usage of the laplace-matrix in example 23
Maxi, citing from the documentation of step-23[1]: "In a very similar vein, we are also too lazy to write the code to assemble mass and Laplace matrices, although it would have only taken copying the relevant code from any number of previous tutorial programs. Rather, we want to focus on the things that are truly new to this program and therefore use the MatrixCreator::create_mass_matrix and MatrixCreator::create_laplace_matrix functions." It comes just handy that we have an easy way to get the required matrices without having to write the assembly loops ourselves. Mutiplying the finite element vector old_solution_u by the laplace matrix is the same as assembling the right-hand side (\nabla phi, \nabla old_solution_u). The point in this example is that you don't have to assemble the laplace matrix for each time step anew. How you do this is totally up to you. In create MatrixCreator::create_laplace_matrix we apply some more optimizations so this might be faster then writing the assembly loop yourself. In the end, you should have a lot of time steps and assembling this matrix once should negligible in comparison to the time you spend in the solver. Hence, my advice would be to use whatever is easier for you. Best, Daniel [1] https://www.dealii.org/8.5.0/doxygen/deal.II/step_23.html#Includefiles Am Freitag, 18. August 2017 17:26:37 UTC+2 schrieb Maxi Miller: > > I am trying to develop my own code based on example 23 (\partial_t U = > i/(2k)\nabla^2U), and am trying to understand the purpose of the > laplace-matrix. Why is it used in order to multiply in both sides, in > comparison with the weak formulation (\nabla\phi_i,\nabla\phi_j) and the > assembly of the usual cell matrix? Assumed I have to split my U into a real > and an imaginary value, is the approach with the Laplace-matrix better, or > the split into the weak formulation and the accompanying matrices on both > sides? > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
