smenon / dynamicTopoFvMesh

Parallel Adaptive Simplical Remeshing for OpenFOAM
http://smenon.github.com/dynamicTopoFvMesh/
38 stars 20 forks source link

3D valve demo: inverting cells #7

Open robertsawko opened 10 years ago

robertsawko commented 10 years ago

Hello Sandeep,

we are trying to move our valve work into 3D. We have some initial success with moveDynamicMesh with a valve moving by a relatively large angle. Here we have an animation.

Unfortunately the case breaks when we go back to sixDoF and pimpleDyMFoam. Here's the configuration with some initial pressure field.

When we use InverseMeanRatio, which was the best in our initial tests, the simulation breaks after 3 time steps with

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 304904 free elements of 305027 total elements.
  This mesh had 305027 tetrahedron elements.
  THERE ARE 5 INVERTED ELEMENTS. 
  (Elements invalid at 5 of 304904 sample locations.)

  5 OF 304904 ENTITIES EVALUATED TO AN UNDEFINED VALUE FOR Inverse Mean Ratio

     Element Quality Statistics

     minimum     average         rms     maximum    std.dev.
-8.808132813e-05 1.173709173 1.185419876 4.7570442710.1662138962

     Number of statistics = 304904
     Metric = Inverse Mean Ratio
     Element Quality not based on sample points.

 Warning: Minimum cell quality is: -0.2921023353 at cell: 298705

 ~~~ Mesh Quality Statistics ~~~ 
 Min: -0.2921023353
 Max: 0.9996758608
 Mean: 0.8661516288
 Cells: 305027

When we use AspectRatioGamma the simulation gets through the first 3 time steps but crashes due to increasing Co number.

Could you please look at the case and advise? I still hope this is some sort of stability problem as for the case that runs we can see pressure oscillating when we admit valve motion. It's just a bit surprising as moveDynamicMesh and pimpleFoam work fine.

smenon commented 10 years ago

This smells like a stability problem related to sixDOF. I would recommend starting with pure mesh-motion (no topology changes) and see how it proceeds. With the right parameters, I'm willing to bet you can go pretty far with that. Once you have that working correctly, you can switch on re-meshing.

robertsawko commented 10 years ago

Thanks for coming back on this issue. Yes, I had a same thought about stability but it seems I didn't really try hard enough. It seems that with slightly ligher valve and accelerationRelaxation parameter set to 0.3 I am getting stable behaviour for AspectRatioGamma. For InverseMeanRatio I need to set the accelerationRelaxation to even lower values, but it's definitely 6DoF which is causing it.

I think this can be closed - we will figure out optimal relaxation values ourselves.

robertsawko commented 10 years ago

I am sorry Sandeep I need to reopen this one. It's clearly a problem with usePointDisplacement. Please look at these two videos

usePointDisplacement setup fixedValuePatches setup

Both are using very similar setups to the one I have on Dropbox already. The calculation done with pimpleDyMFoam and the flow is developing properly. The problem is that when we run in 3D with usePointDisplacement then the points that are not on the patch are not updated in quite the same way.

You can see that fixedValuePatches works perfectly although it keeps the angular momentum constant since in this case it rereads it from a file. The usePointDisplacement version changes the internal points very little and eventually the cells near the moving boundary get squashed. This is probably what was causing negative volume issues.

robertsawko commented 10 years ago

I went back to my rewriting hack which I described in issue #4 i.e. use fixedValuePatches but instead of reading from the dynamicMeshDict I have a separate file that I read and write on each time step that keeps the state of 6DoF. You can see the result here and it seems that the mesh motion works nicely together with pimpleDyMFoam. The cells move all together and the lengthscales of cells adjacent to the valve are more or less constant.

I believe that this shows that there's still something wrong with usePointDisplacement. I am surprised though that it works for 2D nicely. It is a different problem than #4 as back then the whole solver was crashing due to mapping. Here it looks as if the algorithm didn't update point position for internal points.

The position doesn't stabilise but this is probably because our solver setup for this isn't great. We're still in a testing phase.

robertsawko commented 10 years ago

After a week of testing, I must sadly report that the hack isn't really working and there are actually two distinct problems

  1. moveDynamicMesh breaks with usePointDisplacement typically after several time steps. It reports negative quality cells.
  2. Pressure instability in the initial steps of the solvers.

To summarise:

robertsawko commented 10 years ago

I have two examples and a picture. Everything is refereing to moveDynamicMesh which means to me that the issue is with the remeshing and more precisely how pointDisplacment is handled. The case that worked is using fixedValuePatches whereas the case that breaks has usePointDisplacement. The patches are the same and the part which is moving is of uncoupledSixDoFRigidBodyDisplacement with minimal setup - just the angularMomentum is non-zero vector.

usePointDisplacement breaks for me at 0.051 with one inverted cell. I believe I visualised this cell in paraview. highly_skewed_cell When I compare the log files the quality criterium values are not quite the same but even up to the break point they agree up to second significant figure. Then, suddenly usePointDisplacement performs a MESQUITE minimization which gives a slightly lower quality min, performs a number of swaps and breaks in the next time step.

I tried a number of different configs for the remaining pointDisplacement boundaries. In particular I tried fixedValue, zeroGradient and calculated.

I am also attaching the cases if you have the time to look at them. My log files are also there. Note that I am using a lot of Info messages as I am study the execution sequence. https://dl.dropboxusercontent.com/u/13129740/cases.tar.gz

robertsawko commented 10 years ago

I hope you don't mind me coming back here every now and then reporting on the progress. I am not working on this project full time, but rather return to it in spare moments.

So today I spent a happy day debugging the issue. My objective was simple see what usePointDisplacement and fixedValuePatches do differenty. I focused on mesquiteMotionSolver class and more precisely on refPoints_ variable. I made small changes to both your code and 6DoF code and moveDynamicMesh now seems to work for both but I am not happy with these changes.

I am not happy mainly with sixDoF changes. There were two and both pretty simple. Please note that when you go with fixedValuePatches the object is constructed on every time step. With usePointDisplacement it is constructed only once. The main difference that I observed is the initialPoints and initialQ (orientation) which in case of fVP is always the current time step points and Q=Identity. In case of uPD it is "what it says on the tin" that is the points from t=0. The changes I made was to make sure that initialPoints are always current points and Q0 (old value of Q) is identity. This way the behaviour is exactly the same as if the object was reconstructed on each time step.

In your code I changed the line in mesquiteSolver from

if (boundaryConditions_valid())
{
    boundaryConditions_().correctBoundaryConditions();
    refPoints_ = basePoints_() + boundaryConditions_().internalField(); 
}

to

if (boundaryConditions_valid())
{
    boundaryConditions_().correctBoundaryConditions();
    refPoints_ += boundaryConditions_().internalField(); 
}

With these changes moveDynamicMesh seems to give very close results in terms of mesh quality criteria. I was running even on debug 5 to see if the algorithm is doing the same thing and I think it does. However, my changes to sixDoF seem critical but awkward at the same time. The moment I come back to reading the old value of Q properly the refPoints_ become very different, very quickly and the case crashes around 50th time step. So I know I am still not quite there. Maybe you will have some insights as to why this is happening.

smenon commented 10 years ago

Robert,

Sorry for being tardy on this, but since you've been trying a lot of things, it's a bit hard for me to follow what the current state of the problem is.

The trouble with your change from:

refPoints_ = basePoints_() + boundaryConditions_().internalField();

to:

refPoints_ += boundaryConditions_().internalField();

is that the displacement-type boundary conditions (of which six-DOF is one) require that displacements be specified as absolute values from an "initial state", as opposed to incremental changes from the previous time-step (which is what velocity-type boundary conditions do). One potential alternative is to re-write the sixDOF boundary conditions to follow a velocity-type approach, in which case your modification above should work.

robertsawko commented 10 years ago

Hi Sandeep,

don't worry about responses. I hope you don't mind me continuously bumping this thread. I am treating it as my diary on this issue. If I don't write it up somewhere I feel I lost my time somehow. Also, I am new to mesh adaptation in CFD and 6DoF so I had to pick up a lot of stuff on the way. Now I have a good understanding of 6DoF implementation in OF and a some undertanding of your code.

Answering your post, yes, I made these changes to sixDoF so that it wasn't calculating pointDisplacement with respect to initial positions so it was all consistent. But this way isn't optimal as you may get slight elongation or shape change due to numerical errors.

Today I made a different fix and came up with a consistent theory of what may have been going wrong with usePointDisplacement. Let me write it up (again) and I would like to know what you think.

  1. I went back to Port-2.3.x branch
  2. Added this line
    if (boundaryConditions_.valid())
    {
        boundaryConditions_().internalField() = refPoints_ - basePoints_();
    }

after MESQUITE smoothening stuff. My simple moveDynamicMesh tests now work well. I only started larger, flow coupled pimpleDyMFoam cases and will see tomorrow if they work. Now let me discuss a bit what I think is going wrong.

Firstly, everything worked for 2D because you don't do smoothening for 2D. solve in mesquiteMotionSolver returns after surface smoothing so MESQUITE stuff isn't even bothered with any optmization task.

Now in 3D you use the smoothening, but when you apply fixed value patches the refPoints are calculated with respect to basePoints and an updated value of pointDisplacement. The problem is that pointDisplacement is updated only on the patches. The internal point values are all zero even though some points may have been displaced from their initial positions through smoothening. Because of this the information on previous smoothening steps is lost and every time MESQUITE has to optimize from the initial state represented by basePoints_. This may be causing these nasty slivers.

The change I introduced simply saves the pointDisplacement due to smoothening so that next time when we do basePoints_+pointDisplacement.internalField() the internal field values start where we left them last time MESQUITE was doing its job.

Does this make any sense? Am I missing something? If this really fixes anything I am happy to create a pull request but it will be rather minor.

smenon commented 10 years ago

Thanks for the explanation. From what you say, it does seem to make sense, although I can't be certain. Now that you have a case running, it would be nice to see what you come up with. Keep me posted on your progress.

robertsawko commented 10 years ago

It seems to be working. Have a look at the video and at the graph of reported angular velocity. It behaves as I would expect i.e. the valve moves a bit until the hydronamicly induced moment is balanced by gravity moment. Also, using paraview integrate variables filter I double checked that the area of the object did not increase.

I will now run some more tests, including edge refinement options and parallel and if there are no side effects I will close this issue. Thanks for solving issue #2 by the way.

Our aim is to do that sort of calculation on the geometry based on a real system. They're called check-valves, I think and we have a CAD model and the mesh for it. We already know that your code works well with moveDynamicMesh.

robertsawko commented 10 years ago

Hi Sandeep,

so I can confim now that the serial and no edge refinement setup works well. I even updated the video on our basic swing valve geometry using advanced paraview animation features. My colleague will do it on the proper geometry when he comes back from holidays.

Meanwhile though I am struggling with parallel and edge refinement. I made a fork finally as I will need to test various things myself. Parallel run of moveDynamicMesh didn't work even before I introduced pointDisplacement modification after smooth step. The error is

[2] 
[2]  Warning: Minimum cell quality is: -0.2192494193 at cell: 42420
[0] 
[0]  Warning: Minimum cell quality is: -0.1143843639 at cell: 43203
[3] 
[3]  Warning: Minimum cell quality is: -0.299958827 at cell: 5208
[1] 
[1]  Warning: Minimum cell quality is: -0.1329260904 at cell: 2010

 ~~~ Mesh Quality Statistics ~~~ 
 Min: -0.299958827
 Max: 0.9991497622
 Mean: 0.935317963
 Cells: 172545

[2]  Face: 87648 :: 3(622 8246 621) Patch: procBoundary2to3 Proc: 2
[2]  Face: 87649 :: 3(622 121 8246) Patch: procBoundary2to3 Proc: 2
[2]  Face: 82943 :: 3(4527 8246 622) Patch: Internal Proc: 2
[2]  >> Edge: 757::(40 115) mapped: (622 8246)
[2]  Face: 641 :: 3(115 42 40) Patch: Internal Proc: 3
[2]  Face: 786 :: 3(40 115 262) Patch: Internal Proc: 3
[2]  Face: 1066 :: 3(40 115 66) Patch: procBoundary3to2 Proc: 3
[2]  Face: 1092 :: 3(40 96 115) Patch: procBoundary3to2 Proc: 3
[0]  Face: 87389 :: 3(8468 523 522) Patch: FRONTWALL Proc: 0
[0]  Face: 89340 :: 3(522 7181 8468) Patch: procBoundary0to2 Proc: 0
[2]  * * * Error in fillTables * * * 
[0]  Face: 84545 :: 3(522 8468 4239) Patch: Internal Proc: 0
[2]  Edge: 53173 :: (622 8246)
[2]  minQuality: -0.1961357373
[2]  Closed: true
[0]  >> Edge: 4112::(963 819) mapped: (8468 522)
[0]  Face: 5581 :: 3(963 819 1063) Patch: Internal Proc: 2
[0]  Face: 8712 :: 3(1167 963 819) Patch: FRONTWALL Proc: 2
[0]  Face: 9785 :: 3(819 963 818) Patch: procBoundary2to0 Proc: 2
[0]  * * * Error in fillTables * * * 
[0]  Edge: 51220 :: (8468 522)
[0]  minQuality: -0.07954530851
[0]  Closed: false
[1]  Face: 87622 :: 3(8247 7094 7168) Patch: procBoundary1to3 Proc: 1
[1]  Face: 87624 :: 3(8247 7168 1088) Patch: procBoundary1to3 Proc: 1
[1]  >> Edge: 5922::(950 853) mapped: (8247 7168)
[1]  Face: 2070 :: 3(950 854 853) Patch: Internal Proc: 3
[1]  Face: 2109 :: 3(853 1256 950) Patch: Internal Proc: 3
[1]  Face: 2110 :: 3(853 843 950) Patch: Internal Proc: 3
[1]  Face: 8197 :: 3(950 853 855) Patch: procBoundary3to1 Proc: 3
[1]  Face: 8205 :: 3(950 856 853) Patch: procBoundary3to1 Proc: 3
[1]  * * * Error in fillTables * * * 
[1]  Edge: 53153 :: (8247 7168)
[1]  minQuality: -0.02983104867
[1]  Closed: true

Again, the issue is pointDispalcement as the same case runs on fixedValuePatches with the same uncoupled SixDoF setup. It's the smoothening again, isn't it? I probably need to compare again what you've written in fixedValuePatches and see if that can be somehow reapplied to usePointDisplacement.

Edge refinement breaks after taking one whole day of refining a single time step. Unfortunately I lost the log to tell you what the error is exactly. I am in progress of reproducing it now.

So all in all and to put things positive persepctive, there's progress. We could actually carry on this way although we sacrifice speed (no parallel) and a bit of accuracy (we cannot refine).

UlfJan commented 9 years ago

Dear Sandepp, dear Rob,

I am studying mechanical engineering and in order to simulate a not known a-priori 6Dof case, I started working with the dynamicTopoFvMesh class.

I thought that it is a good idea at first to stick to the main branch. So I started with FOAM-extend 3.1 and a usePointDisplacement patch I have found on the Internet(the link to the file is https://www.dropbox.com/s/rbt9lcuq0cjmy50/mesquite.diff.patch?dl=0) .

Thanks to your cases, I managed to create a small test case(the link to the file is https://www.dropbox.com/s/nx8ue85yt2r3uzy/cylinder.tar.gz?dl=0). By running it with moveDynamicMesh and uncoupledSixDoFRigidBodyDisplacement, I think I experienced the exact problem you mentioned before. As the center of mass in pointDisplacement is correct, but by checking in paraview one can see that the position of the cylinder does not match.

"... The trouble with your change from: refPoints = basePoints() + boundaryConditions().internalField(); to: refPoints += boundaryConditions_().internalField(); ..."

And in the patch I used it states exactly this line. That’s why I´m thinking that maybe I do not have the latest version of the usePointDisplacement patch.

My question is, if the usePointDisplacement approach is known to work with foam-extend? Or would it be better to change to another version to continue working on. And if you could please help me out with the usePointDisplacement patch. I would really like to contribute as best as I can. Any help or guidance would be highly appreciated.

Kind regards

robertsawko commented 9 years ago

Hello Jan,

I didn't have the time to work on that issue since I posted here, but using my fork I was able to run several serial cases. I still believe that the explanation I gave above is correct. See my fork as the correction is really simple. Just three lines in an appropriate place.

The only problem is that it breaks parallel which is an issue as our cases take a lot to compute now. I am interested in fixing it, I haven't been updating the code for a while so perhaps with latest fixes it will work out of the box. Within the next three weeks I should actually have some spare time to get back to the parallel issue so maybe we coul take that further.

As for your question you could try running it in foam-extend. For me is not an option as our local HPC runs only canonical OpenFOAM, but if you could confirm that it works there perhaps we could play "spot the difference" game.

smenon commented 9 years ago

Robert and Jan, I've been working on this issue a fair bit recently, but I'm still not close to resolving it, unfortunately. A lot of the recent commits in the Port-2.3.x branch are related to the error that you're seeing here. I know what the problem is - the pointDisplacement pointVectorField is being mapped incorrectly in parallel (and I'm not sure how it's working in serial either). The mapping issue is due to buggy point-mapping in OpenFOAM for pointPatchFields - it's incorrect, mainly because the code there is incomplete. I've been trying to work around the issue by introducing the mesquiteMapper class, but I'm not quite there yet.

I also haven't ported any of this stuff to the extend / master branch either, so I guess it'll have to be with OpenFOAM-2.3.x for now

UlfJan commented 9 years ago

Sandeep and Robert, thank you both very much for your prompt response. After switching to the OpenFOAM-2.3.x branch my test case works with pointDisplacement and moveDynamicMesh. So I have decided to continue my work with OpenFOAM-2.3.x. In order to gain a better understanding, I am going to start by creating a small case for validation using one core. I will keep you updated as soon as there are any news.

smenon commented 9 years ago

Hello Robert and Jan, Could either of you pull the latest fixes from the Port-2.3.x branch and re-test your cases? I believe I've resolved most of the issues that you've been seeing here, except some outstanding stuff related to parallel load-balancing. But point field mapping / edge-refinement / etc should work just fine in parallel. Thanks, Sandeep

robertsawko commented 9 years ago

Great news - thanks. I'll need to recompile the latest OF 2.3.x then as I don't have it currently on none of my machines, but yes we'll test start testing first thing tomorrow morning.

Robert

For those seeking to enter the felicity of matrimonial alliance n/e is the answer. http://en.wikipedia.org/wiki/Secretary_problem

UlfJan commented 9 years ago

Sounds great, thank you very much for the update. I´ll definitely try to run some tests today.

Ulf

robertsawko commented 9 years ago

Ok so I recompiled latest OF 2.3.x and your branch Port-2.3.x. Using scotch decomposition in one my cases ends up with: [3]

[2] 
[2] 
[2] --> FOAM FATAL IO ERROR: 
[2] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[2] 
[2] file: .sectionB_to_sectionA from line 0 to line 0.
[2] 
[2]     From function dictionary::lookupEntry(const word&, bool, bool) const
[2]     in file db/dictionary/dictionary.C at line 442.
[2] 
FOAM parallel run exiting
[2] 
[1] 
[1] 
[1] --> FOAM FATAL IO ERROR: 
[1] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[1] 
[1] file: .sectionB_to_sectionA from line 0 to line 0.
[1] 
[1]     From function dictionary::lookupEntry(const word&, bool, bool) const
[1]     in file db/dictionary/dictionary.C at line 442.
[1] 
FOAM parallel run exiting
[1] 
[0] 
[0] 
[0] --> FOAM FATAL IO ERROR: 
[0] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[0] 
[0] file: .sectionB_to_sectionA from line 0 to line 0.
[0] 
[0]     From function dictionary::lookupEntry(const word&, bool, bool) const
[0]     in file db/dictionary/dictionary.C at line 442.
[0] 
FOAM parallel run exiting
[0] 
[3] 
[3] 
[3] --> FOAM FATAL IO ERROR: 
[3] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[3] 
[3] file: .sectionB_to_sectionA from line 0 to line 0.
[3] 
[3]     From function dictionary::lookupEntry(const word&, bool, bool) const
[3]     in file db/dictionary/dictionary.C at line 442.
[3] 
FOAM parallel run exiting
[3] 
[sow763758c:02152] 3 more processes have sent help message help-mpi-api.txt / mpi-abort
[sow763758c:02152] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages
[2] 
[2] 
[2] --> FOAM FATAL IO ERROR: 
[2] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[2] 
[2] file: .sectionB_to_sectionA from line 0 to line 0.
[2] 
[2]     From function dictionary::lookupEntry(const word&, bool, bool) const
[2]     in file db/dictionary/dictionary.C at line 442.
[2] 
FOAM parallel run exiting
[2] 
[1] 
[1] 
[1] --> FOAM FATAL IO ERROR: 
[1] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[1] 
[1] file: .sectionB_to_sectionA from line 0 to line 0.
[1] 
[1]     From function dictionary::lookupEntry(const word&, bool, bool) const
[1]     in file db/dictionary/dictionary.C at line 442.
[1] 
FOAM parallel run exiting
[1] 
[0] 
[0] 
[0] --> FOAM FATAL IO ERROR: 
[0] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[0] 
[0] file: .sectionB_to_sectionA from line 0 to line 0.
[0] 
[0]     From function dictionary::lookupEntry(const word&, bool, bool) const
[0]     in file db/dictionary/dictionary.C at line 442.
[0] 
FOAM parallel run exiting
[0] 
[3] 
[3] 
[3] --> FOAM FATAL IO ERROR: 
[3] keyword sampleMode is undefined in dictionary ".sectionB_to_sectionA"
[3] 
[3] file: .sectionB_to_sectionA from line 0 to line 0.
[3] 
[3]     From function dictionary::lookupEntry(const word&, bool, bool) const
[3]     in file db/dictionary/dictionary.C at line 442.
[3] 
FOAM parallel run exiting
[3] 
[sow763758c:02152] 3 more processes have sent help message help-mpi-api.txt / mpi-abort
[sow763758c:02152] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

I have usePointDisplacement and a 6DoF type boundary in this case. Can make the file available although would prefer to prepare a minimal example first. I am separately testing edgeRefinement as before it wouldn't work for me even in serial (smoothing error and inverting cells we discussed in #7 ).

kesanyam commented 9 years ago

Hello Sandeep,

I am colleague of Robert and have been using your code for our project.

I have tried out your recent update on branch Port-2.3.x.

The error that Robert had appears to be due to the boundary of the mesh. Because the mesh was created from splitting from a larger mesh using SplitMeshRegions, maybe its configuration on the new boundary is different from the normal boundary. Surprisingly, this isn't an issue for serial run.

I tested your code on a clean mesh on parallel, it runs well for few time-steps.

However, I have a problem when I try to reconstruct the solution. This is the error I got:

--> FOAM FATAL IO ERROR: cannot find file

file:/pimpleDyMFoam_70F/processor0/2.70759e-07/polyMesh/pointProcAddressing at line 0.

From function regIOobject::readStream()
in file db/regIOobject/regIOobjectRead.C at line 73.

FOAM exiting

smenon commented 9 years ago

I think you'll have to use reconstructParMesh first, which matches points geometrically, and then you reconstructPar

kesanyam commented 9 years ago

Hi Sandeep, Thanks for the advice. It works now. The problem with the split mesh can be resolved by altering the boundary file in the polymesh to match those in clean mesh. Yam

kesanyam commented 9 years ago

I have the following error when I attempt to start a new run from the reconstructed solution from pimpleDyMFoam (The simulation crashes and I was trying to run it again). I use adjustableRunTime to save solution at fixed time interval.

--> FOAM FATAL ERROR: Attempt to cast type N4Foam5token8CompoundINS_4ListIdEEEE to type N4Foam5token8CompoundINS_4ListINS_6VectorIdEEEEEE

From function dynamicCast<To>(From&)
in file /home/p004208/OpenFOAM//OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/typeInfo.H at line 93.

FOAM aborting