> On 10/16/2017 09:23 AM, Mark Ma wrote:
> >
> > It almost looks to me like you're applying both Dirichlet values (via
> > the ConstraintMatrix) and Neumann boundary values (via a boundary
> > integral). But then I haven't taken a look at the code, so I can't
> > really say for sure.
> >
> >
> > I think I only applied Dirichlet values via ConstraintMatrix. There
> > is no boundary integrals, so I believe Neumann BC is not
> > implemented.
>
> My guess came from the fact that a Neumann boundary condition would act
> like a boundary heat source. It looks like you have a heat source in the
> last layer of cells around the boundary.
>
> I don't know whether that is true, but it's a place I'd try to
> investigate in debugging.
>
>
Yes, it looks like so. But then I comment the time updating code of rhs
source by the following code, this boundary heat still appears(see figs
below)
//
//time dependent
//assign right hand side
mass_matrix_T.vmult (system_rhs, old_solution_T_cal);
laplace_matrix_T.vmult (tmp, old_solution_T_cal);
system_rhs.add (-time_step * (1-theta), tmp);
/*
assemble_rhs_T (time);
forcing_terms = dynamic_rhs_T;
forcing_terms *= time_step * theta;
assemble_rhs_T (time - time_step);
tmp = dynamic_rhs_T;
forcing_terms.add (time_step*(1-theta),tmp);
system_rhs.add (1,forcing_terms);
*/
//assign system matrix
system_matrix.copy_from (mass_matrix_T);
system_matrix.add (time_step * theta, laplace_matrix_T);
I think this problem lies in the time updating of solution using
old_solution, since the mass_matrix and laplace_matrix have already
eliminated the constraint node,
*mass_matrix_T.vmult (system_rhs, old_solution_T_cal**);* is no longer
valid for this case. This may result into ZEROs at all constrained nodes,
so at boundaries these zero values then propagate into center, as can be
seen from the figure below. For time updating of solution using
old_solution, maybe tutorial-32 could give some ideas. However to
completely understand step-32 may spend some more time for me.
<https://lh3.googleusercontent.com/-yfvZmZDrpHE/WeUMm6TyO-I/AAAAAAAAADo/b0oiTCd5_HMGIsAe9LlZ-RPYF5nWu9ygwCLcBGAs/s1600/tem_no_source.gif>
>
> > But may be you are right since for most of the case, heat equations
> > are solved assuming an adiabatic process, that is dT/dx = 0 at
> > boundaries. This zero values is automatically satisfied.
>
> That's a different question. Whether one may *physically* want a
> different equation is besides the point -- you have chosen one equation
> you'd like to solve, but the numerical solution is wrong. It's
> worthwhile trying to understand why that is before thinking about
> *other* equations.
>
>
> > I think the problem comes from the way of setting Boundary values
> > using constraint method, which in dealii it makes some kind of
> > constraints at vertex of boundaries, i.e. x22 = x1/2 + x2/2,
>
> Hm, that would be the constraint of a hanging node. Do you have hanging
> nodes in your problem? If you do, may I suggest to simplify the problem
> and see if you have the same problem on uniform meshes?
>
>
there is no handing node in my code since I used a uniformed mesh. And also
I tried this, results show no difference.
I am really appreciate your patience to look for this question during your
busy schedule.
Thanks and have a nice day,
Mark
-------------------------------------------------
Laboratoire Hubert Curien, UMR CNRS 5516,
Bâtiment F 18 Rue du Professeur Benoît Lauras
42000 Saint-Etienne
FRANCE
--
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.
