Dear Wolfgang,
Sorry to raise this topic again. But when I add the boundary term to the to
cell_rhs on assembly_system, I receive the following message, based on
which most likely this addition leads to an invertible matrix:
An error occurred in line <457> of file
</usr/local/include/deal.II/lac/solver_cg.h> in function
void dealii::SolverCG<VectorType>::solve(const MatrixType&,
VectorType&, const VectorType&, const PreconditionerType&) [with MatrixType
= dealii::SparseMatrix<double>; PreconditionerType =
dealii::PreconditionIdentity; VectorType = dealii::Vector<double>]
The violated condition was:
false
Additional information:
Iterative method reported convergence failure in step 25. The residual in
the last step was 1.55891e-08.
This error message can indicate that you have simply not allowed a
sufficiently large number of iterations for your iterative solver to
converge. This often happens when you increase the size of your problem. In
such cases, the last residual will likely still be very small, and you can
make the error go away by increasing the allowed number of iterations when
setting up the SolverControl object that determines the maximal number of
iterations you allow.
The other situation where this error may occur is when your matrix is not
invertible (e.g., your matrix has a null-space), or if you try to apply the
wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric
or not positive definite). In these cases, the residual in the last
iteration is likely going to be large.
This is while if I change the pressure, naming it as U2 to another scalar
field variable , naming it as U1, it works, though I get zero iteration for
U1, and in fact it does not solve for U1.
For example if I write:
fe_face_values[U1].value(i, q_point) *... (U1: second scalar variable)
It leads to the error, but if I write:
fe_face_values[U1].value(i, q_point) *... (U1: first scalar variable)
does a number of iterations for U1 , and 0 iteration (probably also
meaning no solution) for U2.
As a bit of further information, my weak form simply ends up to:
dt(M+A)*U1 + A'* U2 = dt*(F_{U1^n-1}) + F_U2
U1 and U2 as two scalar unknowns.
And, the boundary term in matrix-vector form supposed to be my F_U2.
I considered the following block matrices and vectors in forming the linear
system (and solver):
(M+A): system_matrix.block(0,0)
A': system_matrix.block(1,1)
U1: solution.block(0)
U2: solution.block(1)
F_U1: system_rhs.block(0)
F_U2: system_rhs.block(1)
F_U1 would be integration over the domain with phi_i_U1* U1^(n_1)
F_U2 would be integration over the boundary with phi_i_U2*neumann_value
cell_rhs(i) += time_step* phi_i_U1 * old_solution_values *
fe_values.JxW(q_point) ///as (F_U1)
and there after boundary term is added as:
cell_rhs(i) += neumann_value *fe_face_values[U2].value(i, q_point) *
fe_values.JxW(q_point) /// as (F_U2)
and the solve function which could be one of the problems is as the
following:
{
tmp += system_rhs.block(0);
SolverControl solver_control(solution.block
(0).size(), 1e-8 * system_rhs.l2_norm());
SolverCG<Vector<double>> cg (solver_control);
cg.solve(system_matrix.block(0, 0), solution.block (0),
tmp, PreconditionIdentity());
std::cout << " " << solver_control.last_step() << " CG
iterations for U1." << std::endl;
}
{
tmp2 += system_rhs.block(1);
SolverControl solver_control2(solution.block
(1).size(), 1e-8 * system_rhs.l2_norm());
SolverCG<Vector<double>> cg (solver_control2);
cg.solve(system_matrix.block(1, 1), solution.block (1),
tmp2, PreconditionIdentity());
std::cout << " " << solver_control2.last_step() << "
CG iterations for U2." << std::endl;
}
I did quite a number of changes in the code to see the outcome but it seems
I got stuck on an apparently obvious problem and unfortunately haven't been
able to figure it out yet.
I am also not sure if simply addition of boundary term to cell_rhs leads to
the correct recognition as system_rhs.block(1) in the solve?
Would you or anyone else please help me with this problem?
Thanks very much,
Ali
On Monday, January 31, 2022 at 8:08:18 PM UTC+1 Wolfgang Bangerth wrote:
> On 1/31/22 11:41, Ali Seddiq wrote:
> >
> > Thanks for a quick reply, and sorry for ambiguity. To be more precise
> > (trying to be) I should say the lines I modified with the aim of
> > creating a so-called "view" do not do that for my intended variable. My
> > question mainly was that if using the [ ] operator (that I added) does
> > the job of identifying the desired field on the faces? Is the usage
> correct?
>
> Yes. The call in your example code considers the first vector component
> of the i'th shape function at quadrature point q_point.
>
>
> > I am also a bit confused as fe.system_to_component_index
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassFiniteElement.html%23a86644fe67824373cd51e9ff7fca94f8c&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C66a9d6d095e7444d971c08d9e4e94b94%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637792512914352482%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=%2BiFE%2FIimkawMXkp54Z7Kci%2BA7AYn%2B6EJg%2BFeTE5LVYQ%3D&reserved=0>(i)
>
>
> > seems to do the same but not for the faces I guess. So not sure which
> > approach should be taken.
>
> system_to_component index is a property of the finite element, not of
> the FEValues or FEFaceValues class. As such, it can be used as well.
> step-8 shows how it is used for FEValues, but the same could be done for
> FEFaceValues.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> 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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/2f0b7f22-3925-43ac-b56e-b2297072c829n%40googlegroups.com.
