CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Mesquite - Adaptive mesh refinement / coarsening? (https://www.cfd-online.com/Forums/openfoam-solving/110813-mesquite-adaptive-mesh-refinement-coarsening.html)

philippose December 21, 2012 07:19

Mesquite - Adaptive mesh refinement / coarsening?
 
Hello and a good day to everyone :-)!

Its amazing how quickly the year 2012 has gone by, and going by the fact that I am actually typing out a post on the forum, it looks like 21 Dec 2012 is in-fact not the "end of the world" ;-)!

I have a quick query regarding the Mesquite implementation in OpenFOAM by Sandeep and team....

Would it be possible to use Mesquite in order to perform a simulation residual based automatic refinement / coarsening of the mesh during a simulation, so that the mesh is automatically adapted as the simulation progresses?

The documentation of the Mesquite Library from Sandia Labs is not very clear about this, but it looks like the library does support residual based adaptation of meshes (http://acts.nersc.gov/events/Worksho...s/Mesquite.pdf - pages 3-5)

Have there been any thoughts or ideas in this direction? As far as I have understood, the current mesh refinement implementation in OpenFOAM (dynamicrefineFvMesh) only supports hexahedral meshes.... or am I mistaken?

Have a wonderful weekend ahead, and of course, a Merry Christmas, and Happy Holidays to the whole OpenFOAM community :-)!

Philippose

deepsterblue December 21, 2012 10:10

Mesquite currently does not support topological changes. The dynamicTopoFvMesh library does. It currently adapts tetrahedral meshes (using edge refinement / collapse) based on a length-scale field, but that can easily be adapted to represent anything - including residuals.

All you would need to do is specify a scalar value at every cell, and a min/max threshold ratio based on which edges are refined.

philippose December 21, 2012 11:21

Hi Sandeep,

Thanks a lot for your response. I was looking into the ITAPS project a little more in-depth, and found that there was a separate dedicated module called "MeshAdapt" for solution based mesh refinement / coarsening.

Coming back to the OpenFOAM implementation.... Wow :-)! I had not looked into the dynamicTopoFvMesh at all for some reason. Thanks a lot for the heads-up. Shall look into this option and get back to you with either results, or more questions :-)!

Is the solution-mapping in dynamicTopoFvMesh also working? In that.... are the variables mapped from the old mesh to the new mesh automatically as part of the dynamicTopoFvMesh library, or does this need to be done separately?

As of now, the one time I actually tried solution based mesh refinement, was done by:
1. Running the simulation in OpenFOAM
2. Using "simpleFoamResidual" to get the residual field.
3. Using this residual field to generate a "Mesh Size File" for Netgen.
4. Meshing the Geometry in Netgen by using the Mesh Size File for the local refinement / coarsening
5. Importing the mesh back into OpenFOAM
6. Mapping the solution from the old mesh to the new one
7. Rerunning the OpenFOAM simulation

I am quite sure this can be done completely within OpenFOAM.... hence the purpose of the post.



Now again switching back to ITAPS.... Have you worked with components of this project other than Mesquite? I was particularly interested in the EBMesh + MeshKit mix, which seems to be capable of creating Cartesian Meshes with local refinement around geometry specifically for use in "immersed boundary method" (or "Embedded Boundary" as ITAPS coins it) simulations.

Would this be relevant for OpenFOAM given the fact that the last Workshop also had a presentation and implementation of the Immersed Boundary Method in OpenFOAM?

Have a nice day!

Philippose

deepsterblue December 21, 2012 11:28

All the sources are currently hosted here:

https://github.com/smenon/dynamicTopoFvMesh

I've implemented conservative solution mapping in the library, so all fields should be mapped automatically without problems. This should work with OpenFOAM-2.1.x as well - take a look at the Port-2.1.x branch.

It's possible that there's an ITAPS implementation for adaptation - I've never looked into it. Ideally, it's best to keep dependencies to a minimum.

philippose December 21, 2012 14:15

Hello Sandeep,

Ah ok. I thought you were referring to the stock version of dynamicTopoFvMesh already present in the OpenFOAM-1.6-ext git repo.

I have just cloned your github repository, and finished compiling it. Everything went through fine, based on the "Install.txt" file in the repo.

Is there anything I need to know before trying to use the library? Would you happen to have a tutorial somewhere which might be useful as a "getting started" setup?

(I know that there is a circCylinder3d tutorial in the OpenFOAM-ext code)

Philippose

deepsterblue December 21, 2012 14:31

The recent OpenFOAM workshop had a tutorial session with a document / presentation by Kyle Mooney. I guess you could look at those.

philippose December 22, 2012 09:18

Hello again Sandeep,

A Good day to you!

So... I have been slowing putting together a modified version of simpleFoam in order to write out a field with the required mesh size derived from the residuals.

I have a couple of questions regarding the dynamicMeshDict settings for the field based refinement. The documentation from the 7th OpenFOAM Workshop by Kyle Mooney does not really touch on this topic.

I looked into the refinement settings used by interDyMFoam, and though the settings in your library are similar, there are some differences and hence, some questions... :-)!

1. What is implied by the setting: fieldLengthScale?
---- Is this something like a "zero" level for the fieldLengthScale at which no mesh changes would happen (neither refinement nor coarsening)?

2. What does the following setting do: meanScale?
---- This could also be interpreted as a "zero" level, or is this used for something else?

I am assuming that the "lower" and "upper" refinement levels are just the min and max for the lengthScale, which would act like saturations on the lower and upper ends for the mesh edge lengths.... right?

Now.... there are also the settings: minLengthScale and maxLengthScale. How do these play a role in the field based refinement?

Also... are the length scales specified in the lengthScale field used directly as the corresponding edge lengths for the remeshing, or is there some kind of scaling / mapping between the field values and the actual edge lengths?

Wishing you a great Saturday ahead!

Philippose

deepsterblue December 25, 2012 23:49

Philippose,

It's been a while, so I had to look at the sources again to see what was going on.

1. For fieldLengthScale, the intention is to specify a mesh length-scale corresponding to the zero-level. This would imply a mesh edge-length at the "alpha" interface in interFoam, for example.

2. For meanScale, the intention is to specify a mesh length-scale that exists at regions sufficiently far away from the zero-level. This would imply a mesh edge-length at the alpha=0 or alpha=1 regions in interFoam.

I suppose these options don't really sound too intuitive, now that you mention it, but I couldn't really think of anything better at the time.

The lower and upper refinement levels are supposed to be equivalent to the one in interDyMFoam, where an integer value is used to define each level of refinement. In this case, there's no strict "unrefinement" operation that recovers the original mesh, but an edge collapse method, which roughly achieves the same effect.

The min/maxLengthScale values correspond directly to mesh edge lengths. It would probably help to take a look at the lengthScaleEstimator class under the dynamicTopoFvMesh directory to see how things are done.

philippose January 3, 2013 09:19

Hello Sandeep,

A Happy New Year to you :-)! Hope you had a great Christmas break.

I apologise for the delay in replying to your post. I must have missed it. Its only today when I tried to find this post again to ask you something more that I saw that you had already replied to my earlier post.

Thanks a lot for the information. I looked into the lengthEstimator code before seeing your post, and I have some thoughts / questions regarding the implementation of the field based refinement.

1. From the code I understand that the field which drives the length scales is a "volScalarField", and hence, the length scale is defined on a per cell basis rather than a per point basis. This is great because this is what I had in mind earlier.

2. In the part where the field is actually applied to the length scales, I think there is a difference in the way you have implemented it compared to what I had in mind...... Basically, you are checking to see if the average field value of the current cell and the neighbour cell is between the specified "lower" and "upper" refinement levels, and if so, the length scale for that cell is set to the "lengthScale" value specified in the dyamicMeshDict.

To see if I have understood it right.... suppose:
... lengthScale = 0.1
... lowerRefinementLevel = 0.01
... upperRefinementLevel = 0.35
... Average field value = 0.2

In this case, the length scale for that cell and its neighbour would be set to 0.1 right?

That would imply, that the actual length scale of the mesh is not proportional to the field... rather... it is either the already existent mesh length scale (i.e. unmodified), or 0.1.

What I had in mind (and I could be wrong... so please correct me if I am) is, that I could specify a field, calculated based on say... the error residual during the simulation, and the lengthScaleEstimator basically proportionally maps the values in the field to the min and max length scales, giving a resulting length scale which mirrors the values in the field.

For example, if:

... minLengthScale = 0.1e-03
... maxLengthScale = 0.25e-03

... lowerRefinementLevel = 0.2
... upperRefinementLevel = 10.0

... Average field value = 2.0

Then, the resulting length scale for that cell (and the neighbour if applicable) should be:
lengthScale = [(2.0 - 0.2)/(10.0 - 0.2)*(0.25e-03 - 0.1e-03)] + 0.1e-03

And field values above and below the upper and lower refinement levels should be saturated to the upper and lower limits.

What do you think? I hope that makes sense?

Do you think something like this would be a good addition to the current implementation of the remeshing library? And if so, how complicated do you think the implementation would be?

Wishing you a great day ahead!

Philippose

deepsterblue January 3, 2013 10:16

Philippose,

What you're saying makes sense. Again, you're free to use any approach that you find convenient. If you look at the dynamicTopoFvMesh class, you have a DynamicList<scalar> member called lengthScale_, which carries a value specified for each cell in the mesh, and this field that is checked during mesh manipulation. In this case, I fill values in them using the lengthEstimator class (check the dynamicTopoFvMesh::calculateLengthScale function for details), but you could easily fill it using any other volScalarField (will need some minor modifications to the function I mentioned).

What I propose is this - if you can find a scaling function that works for you, create a volScalarField at the top-level in your solver, and add that to the objectRegistry using an IOobject. Once that is done, you can add a dictionary entry to dynamicTopoFvMesh which takes the name of your registered field, and fills values for lengthScale_ from its internal field, rather than calling lengthScaleEstimator for it. Of course, the responsibility of maintaining that field at every time-step is yours. Since it is a registered volScalarField, the mesh will automatically resize it when the mesh changes, so you can be guaranteed that it is properly sized.

Savvy?

philippose January 3, 2013 10:55

Hello again Sandeep,

Cool :-)! Thanks a lot for those tips. I was under the impression that the lengthScaleEstimator was a mandatory part of your system, which cannot be replaced as such.

What I have done so far, is to create a modified version of simpleFoam which writes an additional "volScalarField" containing what in future would be an indication of the length scale for the mesh based on the simulation.

I could change that to make it not just an indicator, but the actual required length scale, which could then be passed on as a field and read in by your library.

I need to have a more detailed peek into the dynamicTopoFvMesh class to determine what the flow of information is like.

Will a change like the one you suggested break any other parts of the library? Or is the lengthScale completely encapsulated and as of now modified only via the lenghScaleEstimator?

In case I have to modify the library, what do you think is the best way to do so? Create a local git branch?


In a way I am actually surprised that no one has wanted a similar solution based re-meshing system for tetrahedral meshes within OpenFOAM so far..... Hmmmm....


Am pretty sure I will be back in a short time with more questions :-)!

Have a nice day ahead!

Philippose

deepsterblue January 3, 2013 11:08

The change I proposed should work. Since the lengthScale_ field is dynamic, values are typically copied or averaged when new cells are added, so I'm assuming that's not a problem for you.

If you're interested in modifying the library, create a fork on GitHub, and once you've tested changes and have it working reliably, initiate a Pull Request, and I can merge it. You have several possible branches to work with - either 1.6-ext or 2.1.x. Pick whichever suits you best.

philippose January 6, 2013 08:27

Hello again Sandeep,

Hope you are having a great weekend!

As you might have already seen, I forked your dynamicTopoFvMesh project in Github and have been working on modifying the code.

I basically finished everything, and tried running a test case using an actual flow solver (variant of simpleFoam) rather than just moveDynamicMesh.

Looks like there might be some issue. In both, the original code, as well as my modified version, the simulation crashes with a segmentation fault when a remeshing is performed (I set the remeshing interval to something like 100, so though the mesh quality and length scales are calculated in each time-step, the actual remeshing only happens every 100 iterations).

After poking around a bit with status outputs at various points in the code, I have found, that the code crashes at the following location:

dynamicTopoFvMesh/fieldMapping/conservativeMapFields.H

between lines 69 and 88.

Surprisingly, it seems to reaching this location whether I have mapping enabled or not (via the "skipMapping" option). I assume that the actual case of whether a mapping is happening or not takes place one level deeper, using either a dummy mapper, or an actual mapping routine.

Just to make it a little clearer.... both the code versions are working fine when the simulation is run with "moveDynamicMesh". It only crashes when executed with an actual fluid solver.

And... I have not made any changes to this part of the code. However, I guess that is more or less clear because the same problem occurs in the original code.


Would it be possible for you to take a quick look into the original code? Or am I missing something?

(P.S.... The actual modification of the length scales seem to be working quite well.... I am still calling "calculateLengthScale". Once things are working you can see whether the approach I have used is acceptable to you or not)

Have a great Sunday!

Philippose

deepsterblue January 6, 2013 11:05

I can take a look. Post a test case at a location that I can access.

philippose January 6, 2013 11:57

Hello Sandeep,

I have copied the test case and the modified version of simpleFoam into my google drive. You can access them using the following link:

https://docs.google.com/folder/d/0B0syP1Z1w2s1YzZiY2UzMTQtNDBhOC00MmRmLWIzZWYtYjJjM TVlYjExZTk3/edit

In the folder, the two relevant files are:

1. mySimpleFoam.tgz - Modified version of simpleFoam to handle dynamic meshes, as well as writing out a field which indicates the length scale

2. 001_Benchmark_Blende_v03_Remeshing_001.tgz - This is a simple test case of an orifice with an inlet and outlet.

In order to run the case, you need to compile "mySimpleFoam", and use this solver.

You can use the original version of the dynamicTopoFvMesh library from your Guthub account without any modifications.

As it is set up now, the remeshing happens in the 10th iteration, and thats when the simulation quits with a segmentation fault, irrespective of whether mapping is enabled or not.


Hope this helps. In case you have trouble accessing the shared files, please let me know.

Philippose

deepsterblue January 6, 2013 12:53

Quick question - did you try this with the Port-2.1.x branch? If not, could you do that and tell me if you get the same error? Unfortunately, I updated my OpenSUSE machine to 12.2, and 1.6-ext no longer compiles with the stock gcc-4.7.1 version, but 2.1.x already has all the source changes in place.

philippose January 6, 2013 13:18

Hello Sandeep,

Nope, I have never used the official OpenFOAM releases yet. Always stuck to the OpenFOAM-ext versions.

I dont know if you have the time / inclination..... I had made a branch on the OpenFOAM-ext Git repository called "hotfix/Gcc47" which I recently updated (merged from the Master) to allow OpenFOAM-1.6-ext to be compiled with GCC-4.7.x (works on my Fedora 17 64-bit VM, which I am currently using).

However, in the meantime I can try to test the library and the case on a version of OpenFOAM-2.1.x once I have the base OpenFOAM-2.1.x version up and running on my system.

Did you try the case on OF-2.1.x with the "mySimpleFoam" solver from me? I dont know if it directly compiles on a 2.1.x version without any changes to the source.

Technically speaking, you can also try the case using "pimpleDyMFoam" in 2.1.x.... you will need to make the relevant changes to "fvSolution" and "fvSchemes" though.

Philippose

deepsterblue January 6, 2013 13:42

I've actually merged the hotfix/gcc47 into master and am compiling right now... I've also tried compiling with 2.1.x, but the errorEstimation directory doesn't exist, and don't know where its been relocated to.

It would be interesting to know if this has any issues in OpenFOAM-2.1.x, just to eliminate the possible compiler-related issues.

deepsterblue January 6, 2013 14:22

4 Attachment(s)
I actually ran it with with master version just now, and everything worked okay... Here's a few screenshots at 0, 10 and 20 (attached), and relevant bits of the log.

Granted, the mesh quality is not that great, but you can play around with that.

philippose January 6, 2013 14:33

Wow... thanks a lot for the information!

Oh ohhhh :-O! Now thats a little worrying.... So basically you ran "mySimpleFoam" on the case I provided you, without any modifications to the case? And it just worked?

Hmmm.... I put in a lot of "Info" statements, and converged on the location of the crash to be the function call: conservativeMapVolFields

Any ideas what the problem on my end could be?

deepsterblue January 6, 2013 14:48

I think I spoke too soon. I forgot to mention that the library was compiled in Debug, and this worked without problems. Now that I compiled optimized, I see the seg-fault. Hmm... Maybe you could run it through valgrind and see?

philippose January 6, 2013 15:02

In a way.... "whew" :-)!

I can/should run it through Valgrind using the optimised build right?

I just deleted all the libraries and applications created by compiling dynamicTopoFvMesh, cleared my "ccache", and shall compile your master branch again from scratch.

In addition, I just cloned the OpenFOAM-2.1.x repository, but getting OpenFOAM-2.1.x compiled, up and running will take a while I think :-)! I am way to used to the "-ext" stream.... specially the Third-party magic my Martin.

Passing information..... I just pushed the modifications I made for the field based scaling into the fork I made.

deepsterblue January 6, 2013 15:32

Okay - the segfault is because of the line:

<< " faceMap: " << fm[faceI] << nl

in topoPatchMapper.C at line: 378. If you comment it out, you get the actual error message, which relates to missing addressing, instead of a seg-fault. As a work-around for now, you can comment out line 363 in topoPatchMapper.C, which is:

if (isA<processorPolyPatch>(patch_.patch()))

This has the effect of artifically mapping from face 0 in the patch, which may not be a big worry in your case, because the problem is steady state. I'll figure out a fix in the mean time.

Learnt a nice lesson today - you can compile optimized with -O3 AND include debug symbols with -g, and the stack trace will contain line numbers instead of plain function names!

philippose January 6, 2013 16:03

Actually, that seems to be what valgrind reported too.....:

The invalid read size of 4 from the calcAddressing function.

Code:

*** Mapping is being skipped ***
 Mapping time: 0.262998 s
 Reordering time: 1.45873 s
==32423== Invalid read of size 4
==32423==    at 0xC1BF3A7: Foam::topoPatchMapper::calcAddressing() const (in /home/philippose/OpenFOAM/philippose-1.6-ext/lib/linux64GccDPOpt/libdynamicTopoFvMesh.so)
==32423==    by 0xC1C0161: Foam::topoPatchMapper::addressing() const (in /home/philippose/OpenFOAM/philippose-1.6-ext/lib/linux64GccDPOpt/libdynamicTopoFvMesh.so)
==32423==    by 0x461C40: Foam::Field<double>::autoMap(Foam::FieldMapper const&) (in /home/philippose/OpenFOAM/philippose-1.6-ext/applications/bin/linux64GccDPOpt/mySimpleFoam)
==32423==    by 0xC06CAFA: void Foam::conservativeMapVolFields<double>(Foam::topoMapper const&) (in /home/philippose/OpenFOAM/philippose-1.6-ext/lib/linux64GccDPOpt/libdynamicTopoFvMesh.so)
==32423==    by 0xC053C55: Foam::dynamicTopoFvMesh::mapFields(Foam::mapPolyMesh const&) const (in /home/philippose/OpenFOAM/philippose-1.6-ext/lib/linux64GccDPOpt/libdynamicTopoFvMesh.so)
==32423==    by 0xC057089: Foam::dynamicTopoFvMesh::resetMesh() (in /home/philippose/OpenFOAM/philippose-1.6-ext/lib/linux64GccDPOpt/libdynamicTopoFvMesh.so)
==32423==    by 0xC05AE85: Foam::dynamicTopoFvMesh::update() (in /home/philippose/OpenFOAM/philippose-1.6-ext/lib/linux64GccDPOpt/libdynamicTopoFvMesh.so)
==32423==    by 0x423A62: main (in /home/philippose/OpenFOAM/philippose-1.6-ext/applications/bin/linux64GccDPOpt/mySimpleFoam)
==32423==  Address 0x0 is not stack'd, malloc'd or (recently) free'd
==32423==
==32423==
==32423== HEAP SUMMARY:
==32423==    in use at exit: 69,190,681 bytes in 856,812 blocks
==32423==  total heap usage: 6,860,255 allocs, 6,003,443 frees, 32,053,601,945 bytes allocated
==32423==
==32423== LEAK SUMMARY:
==32423==    definitely lost: 0 bytes in 0 blocks
==32423==    indirectly lost: 0 bytes in 0 blocks
==32423==      possibly lost: 30,450,860 bytes in 725,719 blocks
==32423==    still reachable: 38,739,821 bytes in 131,093 blocks
==32423==        suppressed: 0 bytes in 0 blocks
==32423== Rerun with --leak-check=full to see details of leaked memory
==32423==
==32423== For counts of detected and suppressed errors, rerun with: -v
==32423== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 2 from 2)
Segmentation fault (core dumped)

Hey thats cool :-)! I had no idea that you can actually get debug symbols on something compiled with a "-O3" !! Always something new lurking around the corner !!

Thanks a lot for the solution idea.... shall do this and let you know what happened in about 45 minutes or so (dinner time now :-)!).

Philippose

philippose January 6, 2013 17:47

2 Attachment(s)
Hi Sandeep,

I made the modification you suggested, and the system now goes past this point without any errors :-)! Thanks a lot.

Now I have come to the next obstacle on the path..... At the iteration when the remeshing takes place, the Poisson flus corrector is called, and then, the usual solver variables are solved.

At this stage, after the "pcorr" Poisson corrector steps, all the residuals in the Ux, Uy, Uz, p, epsilon, k variables jump up by a huge value, and the simulation aborts when solving the epsilon equation.

I have attached a copy of the last few iterations before this takes place.

I tried different settings for the PISO non-orthogonal correctors in the flux corrector step, but nothing seems to work.

Again.... am I missing something, am I doing something wrong, or have you seen this before?

Is there something special to be considered when using the Poisson corrector? Or is it meant for only some specific cases (only laminar, etc..?)?

Also attached is a copy of the checkMesh at the beginning of the simulation.

Philippose

philippose January 6, 2013 17:50

I also tried different solvers and tolerances for the "pcorr" variable required by the flux corrector.... GAMG, and PCG.

Still no change to the situation.

deepsterblue January 7, 2013 09:46

So the debug version of the library seems to proceed normally after remeshing, while the optimized version fails. Could you compare the two by printing out fields before and after field remapping? That might narrow down the scope of the problem. I sure hope this isn't a compiler optimization issue - that would be really annoying...

philippose January 7, 2013 12:13

Hello Sandeep,

A good day to you!

Sorry for the delay in replying to your post. Work started again today.

Regarding the fact that the simulation seems to work fine in debug mode... When you tested that out, did you also have mapping enabled?? The settings in the case I sent you had "skipMapping" set true. The current issue is with "skipMapping" set to false.

Basically, after the discussion we had yesterday, I made the changes suggested by you and then tried to run the simulation with mapping enabled. The code does not fail anymore at the point it used to fail before, but somewhere after the Poisson corrector.

I am wondering whether I should try running the remesher at a higher frequency. It was set to an interval of 200... Maybe something like 20 would make the difference between the original mesh and the remeshed one lower or more acceptable to the flow solvers.

I cannot really compare fields before and after the remeshing because it aborts at a point where none of the new fields have been written yet. I have to see if I can get it to write out "p" and "U" before it tries to solve the turbulence variables.... that way we can compare part of the fields after remeshing.

shall get back with some results once I am back home.

Philippose

danvica January 7, 2013 14:33

Hello Philippose and Sandeep,
Sorry to interrupt your discussion but will the method Philippose presented

Quote:

As of now, the one time I actually tried solution based mesh refinement, was done by:
1. Running the simulation in OpenFOAM
2. Using "simpleFoamResidual" to get the residual field.
3. Using this residual field to generate a "Mesh Size File" for Netgen.
4. Meshing the Geometry in Netgen by using the Mesh Size File for the local refinement / coarsening
5. Importing the mesh back into OpenFOAM
6. Mapping the solution from the old mesh to the new one
7. Rerunning the OpenFOAM simulation

I am quite sure this can be done completely within OpenFOAM.... hence the purpose of the post.

be a generalization of dynamicTopoFvMesh class ?

In specific: can it be used to work with crushing cells (valve closure) ?

I understood Philippose was speaking about refining of a mesh, but what if you send back to Netgen the topological changed mesh ? Can Netgen manage the valve closure case ?

I'm sorry about my question but I must admit I read your discussion and I'm a bit lost.

Thanks.

philippose January 7, 2013 18:10

Hello Sandeep,

So.... I tried the following:

1. Increasing the frequency of the remeshing routine to every 2 intervals. The simulation still aborts with a segmentation fault when solving the epsilon equation while solving for the usual simpleFoam flow variables after the mapping+flux correction (pcorr)

2. Interestingly, I just tried a run with valgrind, and it turned up with the reason for the crash to be:

Code:

==22525== Stack overflow in thread 1: can't grow stack to 0x7fe801f38
==22525== Can't extend stack to 0x7fe800ff0 during signal delivery for thread 1:
==22525==  no stack segment
==22525==
==22525== Process terminating with default action of signal 11 (SIGSEGV)
==22525==  Access not within mapped region at address 0x7FE800FF0
==22525==    at 0x4DAF8C3: Foam::incompressible::RASModels::epsilonWallFunctionFvPatchScalarField::updateCoeffs() (in /opt/Simulation/OpenFOAM/OpenFOAM-1.6-ext/lib/linux64GccDPOpt/libincompressibleRASModels.so)
==22525==  If you believe this happened as a result of a stack
==22525==  overflow in your program's main thread (unlikely but
==22525==  possible), you can try to increase the size of the
==22525==  main thread stack using the --main-stacksize= flag.
==22525==  The main thread stack size used in this run was 8388608.
==22525== Stack overflow in thread 1: can't grow stack to 0x7fe801eb1
==22525==
==22525== Process terminating with default action of signal 11 (SIGSEGV)
==22525==  Access not within mapped region at address 0x7FE801EB1
==22525==    at 0x4801690: _vgnU_freeres (vg_preloaded.c:58)
==22525==  If you believe this happened as a result of a stack
==22525==  overflow in your program's main thread (unlikely but
==22525==  possible), you can try to increase the size of the
==22525==  main thread stack using the --main-stacksize= flag.
==22525==  The main thread stack size used in this run was 8388608.
==22525==

Looks like the stack got filled up..... could it be that the memory is not being freed up? Or is this another one of those misleading errors which has some other root problem?

Have you tried this case on your machine with an optimised build, and with mapping enabled?

Have a great evening ahead, and thanks a lot for the help :-)!

Philippose

philippose January 7, 2013 18:55

I just tried with a "laminar" turbulence model, so that no additional variables other than "p" and "U" are solved for, and now the simulation seems to be going through the remeshing+mapping stages without any problems.

The first run I tried was with the original settings of 100 l/min, and it crashed after 150 iterations, but that may have been because the flow was not laminar.

Its now running with a flow of 50 l/min, and is currently at around 1200 iterations, with remeshing being done every 200 iterations.....

The convergence is not very good though.....

I think some images of the results would be ready tomorrow.

Once unrelated question though..... when you sent me the screen-shots from Paraview the other day, you had somehow managed to cleanly cut the mesh, with the edges only of the cut surface visible instead of the usually messy "surface with edges" result one sees when looking at a "clipped" mesh.... How did you do that??

Have a great evening ahead!

Philippose

deepsterblue January 7, 2013 23:03

Bad convergence could be due to mesh quality. You should have the mesh smoother do more work at each iteration to improve the mesh. Mesh quality typically degrades with a large number of bisections / collapses - you can limit that using maxModifications and / or the sliverThreshold (0.35 is what I typically use). You appear to have a minimum quality of 0.1 after remeshing - which is poor.

Note that you don't "have" to use conservative remapping - take a look at the dynamicTopoFvMesh::mapFields member - it conservatively remaps only for scalars / vectors, but not tensors. You can revert to it to the regular non-conservative way and see if that changes something.

The clipping bit is easy - use the "select cells" widget on the ParaView toolbar after aligning your view to one of the standard axes (can't remember which). Carefully select through using the frustum, invert the selection using the Selection Inspector, and use Filter -> Extract Selection...

philippose January 8, 2013 18:21

Hello Sandeep,

A Good day to you :-)!

I have been looking around through the library today, and have a couple of points.... might be best done in a numbered fashion:

1. First and foremost, thanks a lot for that trick in Paraview for getting a clean cut of the mesh. Things looks so much better and clearer this way :-)!


2. I changed the sliverThreshold factor to 0.35. The mesh smoothing was already set to an interval of "1". I had not looked too deep into the mesh quality outputs during the simulation... sorry about that.


3. Could you give me some more insight into the settings for the smoother? As of now, I have the following primary settings.... can you tell me if they are ok?

========================
tcInner : absGradL2 = 1.0e-04
tcInner : <no limit on CPU time... commented out>

Surface Smoothing CG Solver : tolerance = 1.0e-08
Surface Smoothing CG Solver : nSweeps = 10
Surface Smoothing CG Solver : relaxationFactor = 0.8

surfInterval = 1
========================

Are there any other settings I need to be looking into??


4. I was looking into the field mapping routine in the library, compared to the "mapFields" function in the OpenFOAM fvMesh class.

Is there a specific reason why you only chose to perform a conservative mapping only for vectors and scalars, and not for the other types?

However, just for the sake of clarity.... conservative mapping is primarily required only for transient simulations right? In the case of a steady state simulation, wouldnt the errors which creep in due to a non-conservative mapping get ironed out by the following iterations? Or am I mistaken on this front?

I could remove the conservative mapping, and put back the usual ones from fvMesh.... or is it possible (does it make sense) to call the conservative mapping members also for the tensor types?


5. Following your idea, I compiled my incompressible RAS turbulence models with a "-g" flag, and then ran the simulation with mapping enabled using valgrind, and got the following segmentation fault:

Code:

~~~ Mesh Quality Statistics ~~~
 Min: 0.324886
 Max: 0.99915
 Mean: 0.838364
 Cells: 42648
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

DILUPBiCG:  Solving for Ux, Initial residual = 0.281394, Final residual = 0.000138178, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.988645, Final residual = 0.00056611, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.990267, Final residual = 0.000575066, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.961406, Final residual = 0.00574013, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.655898, Final residual = 0.00285498, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.168068, Final residual = 0.000662986, No Iterations 3
time step continuity errors : sum local = 10.3702, global = 0.00443729, cumulative = 0.00757052
bounding epsilon, min: -14862.6 max: 403126 average: 12662
==27693== Stack overflow in thread 1: can't grow stack to 0x7fe801f58
==27693== Can't extend stack to 0x7fe801010 during signal delivery for thread 1:
==27693==  no stack segment
==27693==
==27693== Process terminating with default action of signal 11 (SIGSEGV)
==27693==  Access not within mapped region at address 0x7FE801010
==27693==    at 0x4DAF8C3: Foam::incompressible::RASModels::epsilonWallFunctionFvPatchScalarField::updateCoeffs() (epsilonWallFunctionFvPatchScalarField.C:173)

The line number 173 in the epsilonWallFunctionFvPatchScalarField.C looks like this:

Code:

void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    // If G field is not present, execute zero gradient evaluation
    // HJ, 20/Mar/2011
    if (!db().foundObject<volScalarField>(GName_)) <-- line 173
    {
        zeroGradientFvPatchScalarField::evaluate();

        return;
    }

    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
    const scalarField& y = rasModel.y()[patch().index()];

Does this ring any bells? It looks like kind of an unlikely place to abort with a segmentation fault unless the object "db()" doesnt exist or "GName_" doesnt exist or something like that right??


On a side note, the minimum mesh quality is 0.324 in this case.... not overly bad I guess given that it started with a minimum quality of 0.308.


So... thats where I have reached so far :-)!

I shall look into running the same simulation without conservative mapping and see if goes through or not....

Wishing you a great evening ahead!

Philippose

philippose January 9, 2013 17:33

Hi Sandeep,

Just wanted to inform you, that I temporarily disabled conservative mapping, and put in the normal non-conservative variant from OpenFOAM for vectors and scalars too.

The simulation now runs even for cases with turbulence enabled. The results are not yet what I am expecting, but there are no more segmentation faults.

Now I need to understand the whole remeshing concept a little more.... and in addition to the questions in the last post, I have a couple more :-)!

1. When I set the "sliverThreshold" to a value which is lesser than the minimum mesh quality as shown by the quality metrics, it causes the remeshing algorithms to be called at every iteration till the minimum mesh quality reaches this threshold, though the "interval" has been set to something else...... Is this behaviour a wanted one?

2. What happens when patches are not specified under any of the types such as "floatingLength", "fixedLength" or "noModification"? How is the lengthScale defined, and under what condition do such patches get refined? The reason for the question is..... it seems that such patches get refined until the "minLengthScale" is reached, whether it is required or not..... is there a specific reason for this behaviour? Or are all patches expected to be categorised into one of the types?

3. In the mesquiteMotionSolver, is the setting "nSweeps" used only for smoothing of patches, or also for smoothing the internal mesh? Is the internal mesh also smoothed each time when the mesquiteMotionSolver function is called from within dynamicTopoFvMesh?


I shall continue looking through the code, and hopefully I will find answers to some of my queries as I go :-)!

Thanks for your help!

Have a great day!

Philippose

deepsterblue January 9, 2013 17:44

1. This is intended behavior. It becomes an issue in transient cases with large deformations, where a large remesh interval is insufficient and the code crashes. This way, a remesh is triggered if the quality goes below the threshold.

2. If no value has been specified, the existing face area is scaled as defined in lengthScaleEstimator::fixedLengthScale. If that is not sufficient, specify a scale that works.

3. The sweeps are only for surfaces. This is basically a linear corrector for curved surfaces. On flat surfaces and the interior mesh, this has no effect. Why would you need sweeps for the interior anyway?

philippose January 9, 2013 18:04

Ah :-)!

1. Ok.... so the sliverThreshold based remeshing is intended.... cool...

2. I shall look into the "fixedLengthScale" member function again, but would you have an idea why these patches (those not explicitly defined) get eventually refined right down to the "minLengthScale"?

3. I was under the impression that the mesh smoothing for the internal mesh and boundary mesh (which I assume is done to improve the mesh quality) is done iteratively, with the mesh quality tending to get better with each sweep, and either stopping when the maximum number of sweeps has been reached, or when say... something like the difference in mesh quality between the previous sweep and the current one goes below some threshold.

It looks like I have probably not understood the way the mesh smoothing has been implemented within the mesquiteMotion library.....


If you have the time, it would be great if you could also look at my post from yesterday? The point 4. in that post is still unclear to me.....

Thanks once again!

Philippose

deepsterblue January 9, 2013 18:17

2. Perhaps you could check the usePolyMesh flag is false (the default) when fixedLengthScale is called, otherwise a scale of 0 is returned, which explains why there's excessive refinement. Be careful though - because it uses the polyMesh face areas, and if this function is called for newly added faces that came from refinement, you'll get a beautiful seg-fault.

3. You may want to set a cpuTime in tcInner options, to ensure that the smoother spends some time smoothing the mesh. This is a mesh optimization process - not a Gauss-Seidel type iterative approach (although that might be happening at the lower levels), so sweeps don't make sense.

No particular reason why it can't be called for tensors - except that you can't take the gradient of a tensor in FOAM - the templates haven't been instantiated yet. One way to enable it is to split tensor fields into 3 vector fields and use a component-wise gradient.

mm.abdollahzadeh January 22, 2013 06:46

Quote:

Originally Posted by deepsterblue (Post 400346)
I actually ran it with with master version just now, and everything worked okay... Here's a few screenshots at 0, 10 and 20 (attached), and relevant bits of the log.

Granted, the mesh quality is not that great, but you can play around with that.

Dear Sandeep

I have used the solver of philippose and its compiled. I have also previously compiled the dynamictopofvmesh. Im working on openfoam 2.1.x . when Im trying to run the test case of philipps, its giving me the following error

Code:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

--> FOAM Warning :
    From function dlOpen(const fileName&, const bool)
    in file POSIX.C at line 1175
    dlopen error : libconservativeMeshToMesh.so: cannot open shared object file: No such file or directory
--> FOAM Warning :
    From function dlLibraryTable::open(const fileName&, const bool)
    in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 96
    could not load "libconservativeMeshToMesh.so"
Create mesh for time = 0

Selecting dynamicFvMesh dynamicTopoFvMesh
Selecting metric Knupp
Selecting motion solver: mesquiteMotionSolver


--> FOAM FATAL ERROR:
solver table is empty

    From function motionSolver::New(const polyMesh& mesh)
    in file motionSolver/motionSolver.C at line 94.

FOAM exiting

Best
Mahdi

deepsterblue January 22, 2013 09:26

Mahdi,

If you want to compile with OpenFOAM-2.1.x, switch to the Port-2.1.x branch of the git repository.

mm.abdollahzadeh January 22, 2013 09:31

Quote:

Originally Posted by deepsterblue (Post 403304)
Mahdi,

If you want to compile with OpenFOAM-2.1.x, switch to the Port-2.1.x branch of the git repository.

Dear Sandeep

I had done that before.
in Fact when i was compiling the dynamictopochange. i had these errors:

Code:

In file included from mesquiteMotionSolver.C:27: mesquiteMotionSolver.H:143: error: ‘Mesquite’ was not declared in this scope mesquiteMotionSolver.H:143: error: template argument 1 is invalid mesquiteMotionSolver.H:155: error: ‘Mesquite’ was not declared in this scope mesquiteMotionSolver.H:155: error: template argument 1 is invalid mesquiteMotionSolver.H:155: error: template argument 1 is invalid mesquiteMotionSolver.H:158: error: ‘Mesquite’ was not declared in this scope mesquiteMotionSolver.H:158: error: template argument 1 is invalid mesquiteMotionSolver.H:161: error: ‘Mesquite’ was not declared in this scope mesquiteMotionSolver.H:161: error: template argument 1 is invalid mesquiteMotionSolver.H:164: error: ‘Mesquite’ has not been declared mesquiteMotionSolver.H:164: error: ISO C++ forbids declaration of ‘TerminationCriterion’ with no type mesquiteMotionSolver.H:164: error: expected ‘;’ before ‘tcInner_’ mesquiteMotionSolver.H:165: error: ‘Mesquite’ has not been declared mesquiteMotionSolver.H:165: error: ISO C++ forbids declaration of ‘TerminationCriterion’ with no type mesquiteMotionSolver.H:165: error: expected ‘;’ before ‘tcOuter_’ mesquiteMotionSolver.C: In member function ‘void Foam::mesquiteMotionSolver::readOptions()’: mesquiteMotionSolver.C:380: error: ‘Mesquite’ has not been declared mesquiteMotionSolver.C:380: error: expected ‘;’ before ‘err’ mesquiteMotionSolver.C:383: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:386: error: ‘Mesquite’ was not declared in this scope mesquiteMotionSolver.C:386: error: template argument 1 is invalid mesquiteMotionSolver.C:388: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:388: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:392: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:395: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:395: error: template argument 1 is invalid mesquiteMotionSolver.C:397: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:397: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:401: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:404: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:404: error: template argument 1 is invalid mesquiteMotionSolver.C:406: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:406: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:410: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:413: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:413: error: template argument 1 is invalid mesquiteMotionSolver.C:415: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:415: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:419: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:422: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:422: error: template argument 1 is invalid mesquiteMotionSolver.C:424: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:424: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:428: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:431: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:431: error: template argument 1 is invalid mesquiteMotionSolver.C:433: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:433: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:437: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:440: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:440: error: template argument 1 is invalid mesquiteMotionSolver.C:442: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:442: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:446: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:449: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:449: error: template argument 1 is invalid mesquiteMotionSolver.C:451: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:451: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:455: error: request for member ‘insert’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:458: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:458: error: template argument 1 is invalid mesquiteMotionSolver.C:460: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:460: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:485: error: request for member ‘found’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:491: error: request for member ‘toc’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:510: error: request for member ‘found’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:514: error: request for member ‘toc’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::qMetricTable_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:613: error: ‘Mesquite’ cannot appear in a constant-expression mesquiteMotionSolver.C:613: error: template argument 1 is invalid mesquiteMotionSolver.C:613: error: template argument 1 is invalid mesquiteMotionSolver.C:613: error: invalid type in declaration before ‘(’ token mesquiteMotionSolver.C:621: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:623: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:623: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:637: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:639: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:639: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:635: warning: unused variable ‘pValue’ mesquiteMotionSolver.C:653: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:655: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:655: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:671: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:673: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:673: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:667: warning: unused variable ‘power’ mesquiteMotionSolver.C:685: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:687: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:687: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:698: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:700: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:700: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:715: error: invalid types ‘int[Foam::label]’ for array subscript mesquiteMotionSolver.C:717: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:717: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:712: warning: unused variable ‘power’ mesquiteMotionSolver.C:745: error: request for member ‘set’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::objFunction_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:747: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:747: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:760: error: request for member ‘set’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::objFunction_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:762: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:762: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:775: error: request for member ‘set’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::objFunction_’, which is of non-class type ‘int’ mesquiteMotionSolver.C:777: error: expected type-specifier before ‘Mesquite’ mesquiteMotionSolver.C:777: error: expected ‘)’ before ‘Mesquite’ mesquiteMotionSolver.C:790: error: request for member ‘set’ in ‘((Foam::mesquiteMotionSolver*)this)->Foam::mesquiteMotionSolver::objFunction_’, which is of non-class type ‘int’


All times are GMT -4. The time now is 01:56.