convergence adaptive refinement

[helmut....@ingenieurbuero-mbi.de]

Hi all,

with the help of you, it was possible for me to complete my code, so I´m able to experiment with different aspects of deal.ii . At the moment I´m capable to solve problems with dirichlet and with robin boundary conditions on meshes build by gmsh read into my program and refined ther (global and adaptive). The program solves the PDEs for thermal conduction and works (partly) as expected, so far, so good.
I examine the convergence of the solution over finer meshes, to get a feeling for the behaviour on different problems.
Now, I have found a behaviour I can not explain or understand.
One problem ( reference case 2 from ISO 10211) does not convergate (right word?) to the correct solution....

I attached a chart with 4 test cases: 2 times a simple square, 2 by 2 m, face 0 and 1 have robin boundary conditions with delta theta 20K. In the corner there is a singularity due to the two different temperatures. In both cases, global and local refinement (kelly) the error gets smaller and smaller, the convergence is as expected, kelly wins in sense of same error with less DOF. The program should be ok.
Btw. as error I assumed the percentage of sum(fluxes) with regard to sum(abs(fluxes) over all boundarys I have set, as I´d assume, that with a perfect solution the sum(flux) should be 0, as I have no sources and the entering flux should be the same as the flux coming out of the other surface.

The other test case (from ISO 10211) behaves quite different: while the global refinement behaves as expected (convergation), the local refinement leads to a growing error. In the uploaded chart there is always a last step with a global refinement of the active cells, so the locally refined ISO-10211 case seems to convergate, but in fact it convergates due to the global refinement in the last step, not due to the local refinement.

I have no idea what happens here, the startig grid is attached too, just for information. (Sorry, but I have no idea how to build a picture of the grid with material IDs). Never the less it is a construction about 0.5m in x-direction, build of 1.5mm steel, wood, 6mm concrete on top and insulation.

Due to the fact that the simple test case behaves as expected, I assume that my code can´t be that wrong. For (other) simple test cases, where I know a analytical solution, the numerical solution ist correct.

Does anybody has an idea what is wrong with the failing test case ?

With best regards,

Helmut

[helmut....@ingenieurbuero-mbi.de]

Ok,

now I have identified at least one problem I have had:

I have done the calculations for integration of fluxes with float variables. The program produces better results with long double, so I assume that there were numerical problems.
Now it is better, but still, below a certain error, the error again rises.

Still under investigation :)

Cheers,

Helmut

[helmut....@ingenieurbuero-mbi.de]

Ok,

now I´m sure that the integration of the flux through the boundary surfaces is correct, as I have build a corrected sum, without further success. I still don´t understand why the solution gets worse during refinement.

Could it be a problem, that the solver builds a solution for my temperature and than, during postprocess the flow is calculated from fe_face_values.grad and the material specific conduction ?

Would it be better to have the flow as unknown, too? Is this possible without to much work ?

Btw, I´m glad, that nobody stated that the solution to my current problem is so simple. As nobody answered at all, I assume that this behaviour is not common.

Cheers, Helmut

[Scott Miller]

Dear Helmut,

Can you be more specific about your problem?  I am not sure what ISO 10211 is; could you point me to a reference on the web that describes the problem?  What equations are you trying to solve?  What discretization are you using, with what finite elements?  These are the sorts of information people need to help you.

Best,

-Scott

[Timo Heister]

> Btw. as error I assumed the percentage of sum(fluxes) with regard to
> sum(abs(fluxes) over all boundarys I have set, as I´d assume, that with a
> perfect solution the sum(flux) should be 0, as I have no sources and the
> entering flux should be the same as the flux coming out of the other
> surface.

Your problems might come from this choice. I would have picked
something like the integral over the square of the flux.


--
Timo Heister
http://www.math.tamu.edu/~heister/

[helmut....@ingenieurbuero-mbi.de]

Hello Scott,

thank you for the reply, you are right, my post was very unspecific, or at least just understandable for people involved in building physics. So I´ll try to explain.

I solve problems in building physics, mainly edges or corners of buildings or other details of buildings. The goal is to decide, if a given construction is faulty (in german : mangelhaft) or not, so in the next step the question can be answered if the user who lives in this building is responsible for moulding or the building itself.

In eurobe there are standards for the necessary calculations, they define how exact the program has to calculate, how exact the results have to be presented and so on. The mentioned DIN EN ISO 10211 gives advice how to cut out the interesting part out of the building and how you build the model for further simulation.

Therefor you have to solve mainly the heat transfer equations for the steady state (constant boundary conditions in time) so i assume the class of problem is called Laplace. The main target ist the minimum surface temperature on the inner boundary with given boundary conditions.

In the DIN EN ISO 10211 are 2 examples that are used to verify computer programms to assure that only validated software is used to assess a given construction.

With the second example i got into trouble. I have made a picture from the model, see attachment. Shown is only the left part of the model. As you can see there a 4 different materials and robin bc on bottom on on top, left and right side are adiabatic.

As a criterion for "good enough" mesh, ISO 10211 defines the fraction of sum(flow throug boundaries) and sum(abs(flow)) has to be less than o given factor.

Now I have played around with adaptive refinement of the mesh, as this was one of the main reasons why I´m interested in deal.ii . While the here shown model is very small ( about 5cm height, 50cm wide) in our work we have often 3-d models 3m by 1.5m by 1.5m or, with earth contact larger models. So it would be very helpful just to start with a coarse grid and have adaptive refinement do the work.

But unfortunately I observed - during playing with the test case from ISO 10211 - that my error increases during refinement, but only with this model while a simple test case - just a square, 1 material, 2 bc - bahaves as expected.

I have attached the gmsh grid, for your information. In my program I have first done an anisotropic refinement to assure that the cells have an aspect ration better than 2, but it doesent matter : with or without the initial anisotropic refinement the convergence does not happen as expected ( see my first post).

The adaptive refinement during itterative calculations is done by kelly error estimator, which should be ok for thermal conductivity, while the initial anisotropic refinement is done by hand : I just compare dimensions in x and in y direction and decide to refine in one direction or not.

The program uses simple finite elemnts (Q1), I tried with (bi-)linear elements and with (bi-)quadratic elements. As expected better results with quadratic elements as they better fit to the problem, but in total the same problem with regard to the flow-error.

And just to be complete, I calculate the flows throug the boundaries in this way, which I assume is correct:

inline double myStep3::postprocess()
    {
    QGauss<1> face_quadrature_formula(3);
    const unsigned int n_face_q_points = face_quadrature_formula.size();

    FEFaceValues<2> fe_face_values(fe, face_quadrature_formula,
        update_gradients | update_normal_vectors | update_JxW_values);

    long double summe_q_boundary = 0.0;
    long double summe_abs_q_boundary = 0.0;

    std::vector<MyBoundary> theBCs;
    theBCs = myBCSet.getAllBC();

    typedef std::vector<MyBoundary>::iterator BC_Itterator;
    for (BC_Itterator it = theBCs.begin(); it != theBCs.end(); ++it)
    {
    int bcNummer = it->get_Nummer();
    long double summe_q_face = 0.0;

    DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(),
        endc = dof_handler.end();
    for (; cell != endc; ++cell)
        {
        // Lambda der Zelle bestimmen
        long double lambda = myMatSet.getLambdaByNummer(
            static_cast<int>(cell->material_id()));

        // alle Faces der Zelle abarbeiten
        for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell;
            ++face)
        if (cell->face(face)->at_boundary()
            && (cell->face(face)->boundary_indicator() == bcNummer))
            {
            // das Face liegt auf der aktuellen Boundary, also integrieren
            fe_face_values.reinit(cell, face);

            // gradienten der Lösung auf der Oberfläche bestimmen
            std::vector<dealii::Tensor<1, 2>> grads(n_face_q_points);
            fe_face_values.get_function_gradients(solution, grads);

            for (unsigned int q = 0; q < n_face_q_points; ++q)
            {
            // Wärmestromdichte = -k*grad(T)
            double summand = -1.0 * lambda * fe_face_values.JxW(q)
                * (grads[q] * fe_face_values.normal_vector(q));
           
            // naive addition, funktioniert aber mit double
            summe_q_face += summand;
            }
            }
        }
    std::cout.setf(std::ios::fixed, std::ios::floatfield);
    std::cout.precision(12);
    std::cout << "Wärmestrom Boundary Nr " << bcNummer << " ist : "
        << summe_q_face << " W/m" << std::endl;

    summe_q_boundary += summe_q_face;
    summe_abs_q_boundary += abs(summe_q_face);
    }

    std::cout << "Summe Wärmeströme            : " << summe_q_boundary << " W/m"
        << std::endl;
    std::cout << "Summe Absolute Wärmeströme   : " << summe_abs_q_boundary
        << " W/m" << "\n" << std::endl;
    std::cout << "Fehler                       : "
        << abs(summe_q_boundary) / summe_abs_q_boundary * 100.0 << "% \n\n"
        << std::endl;

    konvergenz << dof_handler.n_dofs() << "\t";
    konvergenz << abs(summe_q_boundary) / summe_abs_q_boundary * 100.0 << "\n";

    return summe_abs_q_boundary;
    }

I´m sorry for this quite long post, and thank you for your attention!

With kind regards,

Helmut

[helmut....@ingenieurbuero-mbi.de]

Hello Timo,

thank you for your reply.

I have to choose those quantities (see my previous post, where I have written why). Within other fem code I have never had those problems before. But to be honest, this is the first time i work with adaptive refinement.

I have examined the parts of the sum, and I see nowhere a reason why this choice should have impact on the mentioned problem, at least with double it works. Please, can you explain that any further?

With kind regards,

Helmut

[Bärbel Janssen]

Hi Helmut,

I'm afraid we are still not able to help you with the information you provided in your long email. There are some questions not answered yet:

* Do you use continuous finite elements (FEQ<2>(order))? If you do that, do you use a ConstraintMatrix for the hanging nodes in adaptive refinement?

* Which solver do you use? To which tolerance do you solve your problems?

* In which sense does your quantity not converge? Does this occur on coarse grids or rather on finer grids?

* What is the specific configuration where you have convergence? (FE, tolerance of the solver, adaptive, no adaptive refinement, isotropic, anisotropic refinement)?

Best,

Bärbel

2013/6/19 <helmut....@ingenieurbuero-mbi.de>

--
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+un...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[helmut....@ingenieurbuero-mbi.de]

Hello Bärbel,

thank you for your attention,

I´m sorry that still relevant facts are missing. Unfortunately I´m not sure which facts are necessary for you, but due to your question I can try to answer, thank you.

As finite elements  I use (quite common?) FE_Q<2>, as stated above linear and also (different test) quadratic.

And I use the constraint matrix to handle the hanging nodes, as mentioned in the tutorials.

The solver is SolverCG with tolerance 1e-12, preconditioned with SSOR, initialized to 1.0 (FYI: this is not a decision I could justify, but working throug the tutorials this appeared into my code)

The solution build by the solver does converge as expected, the unknown temperature looks right, so I think this part works like expected. But the derived quantity (heat flow, build by conductivity and the gradient of the solution) integrated over the boundaries does not converge, at least for this example. For the simple square domain this error get smaller and smaller (chart in my first post).

On the uploded domain representing the DIN EN ISO 10211 example as you can see in the chart from the first post the error grows while adaptive refinement increases number of Dof. The error then gets smaller when the locally refined grid gets refined globaly (last plotpoint in chart).

To your last question : all 4 curves shown in the chart from the first post were build with the same code, only an other grid was used and the boundary conditions were adjusted.

I hope this helps to better understand me. If it would help, I could post my complete code, but to compile it a few further helper classes would be required (simple dataobjects for material and boundaries and their collections) I could upload them too, but I don´t want to flood you with the code...

With kind regards,

Helmut

[Bärbel Janssen]

Hi Helmut,

thank you for making it more clear for me. One more question: Did you check that all the desired faces on the boundary are integrated as expected? How do you colorize your grid? Is the colorization inherited in the adaptive refinement?

Best,
Bärbel

[helmut....@ingenieurbuero-mbi.de]

Hi Bärbel,

I´m sorry that you and all other have to ask for information I should have given already...

Yes, at first view it looks like all faces are integrated. Unfortunately there are lot of faces so it is not easy to be sure. I´ll think about a possibility to check it.

I read in a colorized mesh build by gmsh. Within gmsh I set the boundary indicators and the material indicators. Due to the fact that both are recognized and the calculated temperature makes sense, I think this should be ok also for the refined grids.

I use gmsh because geometries I will have to calculate in the future are to complex to be build by hand or with the primitives in deal.ii . In the past I have written different "geometrie-generators" that produce geo-files for gmsh. At the moment I do the calculations with getdp that understands gmsh-meshes. With small modifications it was possible to produce the grids that deal.ii can handle.

With best regards from hot germany,

Helmut

Am Mittwoch, 19. Juni 2013 14:20:03 UTC+2 schrieb Bärbel Janssen:

Hi Helmut,

thank you for making it more clear for me. One more question: Did you check that all the desired faces on the boundary are integrated as expected? How do you colorize your grid? Is the colorization inherited in the adaptive refinement?

Best,
Bärbel

[helmut....@ingenieurbuero-mbi.de]

Hi,

is there a way to output a (refined) grid with colorization (eg. with boundary-id and material-id as number) just to be sure that the inheritance is performed properly ? Which other way would be thinkable ?

Cheers,

Helmut

Am Mittwoch, 19. Juni 2013 14:20:03 UTC+2 schrieb Bärbel Janssen:

[Scott Miller]

is there a way to output a (refined) grid with colorization (eg. with boundary-id and material-id as number) just to be sure that the inheritance is performed properly ? Which other way would be thinkable ?

You can do something like this to output material ids:

Vector<int> material_ids (triangulation.n_active_cells());

// Omitted code:  Create iterators:  cell, endc

unsigned int cell_counter =0;

for (unsigned int cell_counter =0; cell!=endc; ++cell,++cell_counter)

    material_ids(cell_counter) = int(cell->material_id());

// Add to DataOut object

data_out.add_data_vector (material_ids, "Material IDs");

I don't know that you can output boundary ids to visualize easily, at least I typically don't.  You can, however, do the same sort of thing but give cells with certain boundary faces different values in the data vector.

-Scott

    

 

[Bärbel Janssen]

You could use GridOut to write it as a ucd file. If you pass the correct flags to GridOut it colorizes the boundaries and you can visualize it with paraview. That's how I usually do it.

There is one more question: How does the solution on the adaptive mesh look compared to the uniform refined mesh? Is there something wrong with the solution or just with the quantity you compare?

Best,
Bärbel

2013/6/19 Scott Miller <miller....@gmail.com>

 

--

[helmut....@ingenieurbuero-mbi.de]

Hello Bärbel and Scott,

thank you for the response.

I´ll give it a try, both tips seem to be helpfull to me. Unfortunately I have a lot of work at the moment, so I´m not shure when I´ll be able to proceed. Never the less, thank you and I´ll post a notice when I know more.

Cheers,

Helmut

[helmut....@ingenieurbuero-mbi.de]

Hi Bärbel,

could you give me a short advice how to export the grid/mesh to ucd as you suggested in a previous post ?

I tried it this way:

stringstream ss2;
    ss2 << "grid-" << cycle << ".ucd";
    GridOut myGridOut;
    GridOutFlags::Ucd flags;
    flags.write_faces = true;
    flags.write_lines = true;
    myGridOut.set_flags(flags);
    std::ofstream output(ss2.str());
    myGridOut.write_ucd(triangulation, output);

but ParaView doesn´t find the mesh :

ERROR: In /Users/kitware/Dashboards/MyTests/NightlyMaster/ParaViewSuperbuild-Release/paraview/src/paraview/Utilities/VisItBridge/databases/AvtAlgorithms/vtkAvtSTMDFileFormatAlgorithm.cxx, line 122

vtkVisItCEAucdReader (0x10b664120): Unable to find any meshes

Any help would be apreciated :)

Cheers,

Helmut

Am Mittwoch, 19. Juni 2013 15:10:45 UTC+2 schrieb Bärbel Janssen:

[Wolfgang Bangerth]

Helmut,
I think the easiest way would be to post the code you currently have, together
with all inputs. Then also describe how you think the solution should look
like and how it actually does look like.

It may be possible to see right away what is going on.

Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Lam DANG]

Dear All

I use write_ucd with flags.write_faces = true, but the result of .ucd file is lack of information of face id 0.
How I can output the information of face id 0?

The piece of code:
    Triangulation<2> triangulation;
    GridGenerator::hyper_cube (triangulation, 0, 1, true);
    GridOut grid_out_ucd;

    GridOutFlags::Ucd flags;
    flags.write_faces = true;

    grid_out_ucd.set_flags(flags);
    std::ofstream out_ucd ("triangulation.ucd");
    grid_out_ucd.write_ucd (triangulation, out_ucd);

The result of "triangulation.ucd"
4 4 0 0 0
1  0 0 0
2  1 0 0
3  0 1 0
4  1 1 0
1 0 quad    1 2 4 3
2  2  line    1 2
3  1  line    2 4
4  3  line    3 4

The information of line 0 (5 0 line 1 3) is not in output file

Best regards

Lam DANG

[Wolfgang Bangerth]

On 06/06/2017 04:43 AM, Lam DANG wrote:
> Dear All
>
> I use /write_ucd/ with /flags.write_faces = true/, but the result of
> .ucd file is *lack of information of face id 0.*

> How I can output the information of face id 0?
>
> The piece of code:

> / Triangulation<2> triangulation;

> GridGenerator::hyper_cube (triangulation, 0, 1, true);
> GridOut grid_out_ucd;
> GridOutFlags::Ucd flags;
> flags.write_faces = true;
> grid_out_ucd.set_flags(flags);
> std::ofstream out_ucd ("triangulation.ucd");

> grid_out_ucd.write_ucd (triangulation, out_ucd);/
>
> The result of "/triangulation.ucd/"
> /4 4 0 0 0

> 1 0 0 0
> 2 1 0 0
> 3 0 1 0
> 4 1 1 0
> 1 0 quad 1 2 4 3
> 2 2 line 1 2
> 3 1 line 2 4

> 4 3 line 3 4 /
>
> The information of line 0 (/5 0 line 1 3/) is not in output file

For reasons that are no longer entirely clear to me, GridOut documents
this behavior as follows:

struct Ucd
{
[...]

/**
* When writing a mesh, write boundary faces explicitly if their
boundary
* indicator is not the default boundary indicator, which is zero.
This
* is necessary if you later want to re-read the grid and want to
get the
* same boundary indicators for the different parts of the boundary
of the
* triangulation.
*
* It is not necessary if you only want to write the triangulation
to view
* or print it.
*
* Default: @p false.
*/
bool write_faces;

In other words, it *purposefully* omits faces with a zero boundary
indicator. In your case, since you are "coloring" the cube upon mesh
generation, the face that's omitted is exactly the one that has boundary
indicator zero.

In other words, the function really does as is intended. What the
original intent of the rule to omit faces with zero boundary indicator
is not completely clear to me any more, however.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/