[Kellis]

Originally Posted by bigphil

Hi Kellis, Dave,

I suggest running the solver through a debugger (e.g. gdb) in debug mode and finding the lines of code where the solver stops. From there, it should be clear where the error is coming from and then hopefully a solution can be found.

Philip

Philip,


Thank you for the response. After recompiling with the Debug option set, and using mpirunDebug to run the case as follows:

mpirunDebug -np 4 icoFsiElasticNonLinULSolidFoam -parallel

I get the following error returned:

[...]

Create fluid-to-solid and solid-to-fluid interpolators
Check fluid-to-solid and solid-to-fluid interpolators
[0] 
[0] 
[0] --> FOAM FATAL ERROR: 
[0]     incompatible fields Field<vector> f1(2842), Field<vector> f2(2842) and Field<vector> f3(0)
    for operation f1 = f2 + f3
[0] 
[0]     From function checkFields(const UList<Type1>&, const UList<Type2>&, const UList<Type3>&, const char*)
[0]     in file /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/FieldM.H at line 84.
[0] 
FOAM parallel run aborting
[0]

By manipulating the decomposeParDict, I have been able to determine that the vector fields f1 and f2 have the same size as the FSI patch assigned to processor0 (i.e, there are 2,842 faces of the FSI patch assigned to processor0 above), and f3 seems to have the same size as the FSI patch assigned to processor1. I'm not sure what operation the solver is trying to do here; do you have any ideas?


Thank you,
Kellis

[bigphil]

Hi Kellis,

Try the following command to run the solver in parallel through the gdb debugger:

mpirun -np 4 xterm -e "gdb --args icoFsiElasticNonLinULSolidFoam -parallel"

This will open up four terminals, where each processor is running in a different terminal.

Once gdb opens, you can start each process with the "r" command. You should then be able to see which processor stops and you can find out where it stopped using the "bt" (backtrace command). You may need to add a breakpoint at FieldM.H:84 to find out where.

Best,
Philip

[Kellis]

Originally Posted by bigphil

Hi Kellis,

Try the following command to run the solver in parallel through the gdb debugger:

mpirun -np 4 xterm -e "gdb --args icoFsiElasticNonLinULSolidFoam -parallel"

This will open up four terminals, where each processor is running in a different terminal.

Once gdb opens, you can start each process with the "r" command. You should then be able to see which processor stops and you can find out where it stopped using the "bt" (backtrace command). You may need to add a breakpoint at FieldM.H:84 to find out where.

Best,
Philip

Philip,


Thanks for your continued help. Per your instructions, I ran the program in gdb, and added a break at that location. The output at the break was as follows:

Breakpoint 1, Foam::checkFields<Foam::Vector<double>, Foam::Vector<double>, Foam::Vector<double> > (f1=..., f2=..., f3=..., op=0x553885 "f1 = f2 + f3")
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/FieldM.H:84
84            )   << "    incompatible fields"

Running the backtrace produced the following output:

#0  Foam::checkFields<Foam::Vector<double>, Foam::Vector<double>, Foam::Vector<double> > (f1=..., f2=..., f3=..., op=0x553885 "f1 = f2 + f3")
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/FieldM.H:84
#1  0x00000000004d07ab in Foam::add<Foam::Vector<double>, Foam::Vector<double> > (res=..., f1=..., f2=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/FieldFunctions.C:689
#2  0x00000000004b4652 in Foam::operator+<Foam::Vector<double>, Foam::Vector<double> > (f1=..., f2=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/FieldFunctions.C:689
#3  0x00000000005008a5 in Foam::sumOp<Foam::Field<Foam::Vector<double> > >::operator() (this=0x7fffffff47c0, x=..., y=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/ops.H:158
#4  0x00000000004e7b89 in Foam::Pstream::gather<Foam::Field<Foam::Vector<double> >, Foam::sumOp<Foam::Field<Foam::Vector<double> > > > (comms=..., Value=..., bop=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/gatherScatter.C:77
#5  0x00000000004cd5bd in Foam::reduce<Foam::Field<Foam::Vector<double> >, Foam::sumOp<Foam::Field<Foam::Vector<double> > > > (comms=..., Value=..., bop=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/PstreamReduceOps.H:48
#6  0x00000000004b2913 in Foam::reduce<Foam::Field<Foam::Vector<double> >, Foam::sumOp<Foam::Field<Foam::Vector<double> > > > (Value=..., bop=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/PstreamReduceOps.H:63
#7  0x000000000049ded5 in main (argc=2, argv=0x7fffffffae68) at createZoneToZoneInterpolators.H:53

Unfortunately I don't have much experience with C++ so most of this is gibberish to me. Are you able to extract anything useful from this info?


Thanks again,
Kellis

[bigphil]

Thanks Kellis,

From the last line of the backtrace, we can see that the code is failing in the file "createZoneToZoneInterpolators.H" at line 53: this line of code is:

       reduce(fluidZoneFaceCentres, sumOp<vectorField>());

which is a parallel syncing step (reduce) where the vectorField "fluidZoneFaceCentres" is summed across all processors. The size of the "fluidZoneFaceCentres" field should be the same on all processors (it is based on the fluid interface globalFaceZone); however, from the error you receive, we can see that the size of this field is different on each processor.
On proc0 we can see that it has size "2842" whereas on proc1 it has size "0".

So it looks like the faceZone is not a globalFaceZone; did you include the FSI fluid interface zone name in the globalFaceZones list within decomposePar?

Philip

[Kellis]

Originally Posted by bigphil

So it looks like the faceZone is not a globalFaceZone; did you include the FSI fluid interface zone name in the globalFaceZones list within decomposePar?

Philip

Philip,

Wow, you are a magician! Indeed I don't currently have globalFaceZones option set in the decomposeParDict. I have tried this before, but it produces the following error while decomposing the fluid fields (full log attached):

Processor 0: field transfer
--> FOAM Warning : 
    From function dlLibraryTable::open(const fileName& functionLibName)
    in file db/dlLibraryTable/dlLibraryTable.C at line 124
    could not load libgroovyBC.so: cannot open shared object file: No such file or directory
Initializing the GGI interpolator between master/shadow patches: side1/side2



--> FOAM FATAL ERROR: 
size of field = 129349 is not the same as the size of mesh = 126492

    From function DimensionedField<Type, GeoMesh>::DimensionedField(const IOobject& io,const Mesh& mesh, const dimensionSet& dims, const Field<Type>& field)
    in file /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/DimensionedField.C at line 71.

Looking into the fluid/processor0/0 directory, I noticed it didn't write the motionU file, so I'm guessing this is where the error lies. I went ahead and ran decomposePar again with gdb, and got the following backtrace:

#0  0x00007ffff1e4b428 in __GI_raise (sig=sig@entry=6)
    at ../sysdeps/unix/sysv/linux/raise.c:54
#1  0x00007ffff1e4d02a in __GI_abort () at abort.c:89
#2  0x00007ffff352ec94 in Foam::error::abort (
    this=0x7ffff401fa20 <Foam::FatalError>) at lnInclude/error.C:252
#3  0x00000000005a733a in Foam::operator<< <Foam::error> (os=..., m=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/errorManip.H:85
#4  0x000000000062de1e in Foam::DimensionedField<Foam::Vector<double>, Foam::pointMesh>::DimensionedField (this=0x93686b0, io=..., mesh=..., dims=..., 
    field=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/DimensionedField.C:66
#5  0x00000000006065cb in Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh>::GeometricField (this=0x93686b0, io=..., 
    mesh=..., ds=..., iField=..., ptfl=...)
    at /home/kckincaid/foam/foam-extend-3.2/src/foam/lnInclude/GeometricField.C:327
#6  0x00000000005d9800 in Foam::pointFieldDecomposer::decomposeField<Foam::Vector<double> > (this=0x7fffffff35f0, field=...)
    at pointFieldDecomposerDecomposeFields.C:130
#7  0x00000000005c03a4 in Foam::pointFieldDecomposer::decomposeFields<Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh> > (
    this=0x7fffffff35f0, fields=...)
    at pointFieldDecomposerDecomposeFields.C:142
#8  0x00000000005b0f0c in main (argc=2, argv=0x7fffffffb498)
    at decomposePar.C:746

Looking at the decomposePar.C file at line 746, I see this is where the pointVectorFields get written, which reinforces the idea that there's an issue with the motionU field. However, I'm not sure where decomposePar is getting the size of the field and mesh from. I searched my case directory for the two size values it spit out, but they don't seem to be written in any of the files.

I've attached the full decomposePar log file, as well as the dictionaries for both the fluid and solid meshes. Hopefully these are helpful.

Again, thank you for all of your assistance so far!

Thanks,
Kellis

[bigphil]

Hi Kellis,

Yep it looks like it is a problem decomposing the motionU file. To double-check, I suggest you move the motionU file out of the "0" directory and see if decomposePar works in that case.
If it does then we know motionU that is causing the problem. motionU is the mesh motion field for the FE motion solver so I think decomposePar doesn't know what to do with it when there are globalFaceZones. I vaguely remember there being a fix for this in decomposePar (maybe it was a later commit of foam-extend-4.0). Which commit of foam-extend-4.0 are you using? You can check by typing "git log" in the foam-extend-4.0 directory (assuming you downloaded it with git). I am using commit 246a172c9d9ecdbfc8f18b8224467de7f2b7606e dated 21 Sep 2017.

Also, can you share a simple version of your case and I can check if it decomposes with version?

Alternatively, you could also try the FV motion solver in the dynamicMeshDict: it uses pointMotionU instead of motionU.

Philip

[Kellis]

Originally Posted by bigphil

Also, can you share a simple version of your case and I can check if it decomposes with version?

Alternatively, you could also try the FV motion solver in the dynamicMeshDict: it uses pointMotionU instead of motionU.

Philip

Philip,

I will try to get a simplified version of my case put together for you and send it to you via PM. In the mean time, I checked if it will decompose without the motionU or pointMotionU files, and indeed it works perfectly.

I did a little more investigation into the decomposition error. As a reminder, it was throwing the following error while trying to decompose with the motionU files in the 0 directory:

--> FOAM FATAL ERROR: 
size of field = 129349 is not the same as the size of mesh = 126492

I ran the checkMesh utility in the processor0 directory, and found out the following:

Mesh stats
    all points:           129349
    live points:          126492
    all faces:            353287
    live faces:           350472
    internal faces:       322452
    cells:                112154
    boundary patches:     12
    point zones:          0
    face zones:           4
    cell zones:           0

So it appears it is complaining about the discrepancy between the total number of points, and the number of live points. I can't find any documentation on what exactly "live" points are; in the full mesh, the number of live points exactly matches the total number of points.

So I suspect it is probably something with the way the mesh is getting decomposed. I'll send you my case ASAP.

As always, thank you!

-Kellis


EDIT: Also, I forgot to mention I am currently using Extend v3.2, commit 56c5fe3e449d4c1c1c1a5f56c9320c56f65785bc, dated 22 Sept 2016. I tried using FE 4.0, but had issues with the solver hanging while setting up the GGI interpolator.

[bigphil]

FYI, for others:

For foam-extend-3.2, using decomposeParFsi from the extendBazaar fixes this issue.

For foam-extend-4.0 (and later), this issue has already been fixed in decomposePar.

Philip

[wulonglong]

Hello all,

Is there any recommended solver for this FSI problem (in attachment)?

I hope me ask in the right thread...

Many thanks in advance.

-Wu

[bigphil]

Hi Wu,

I have been working on this "flexible" dam break case using the solids4foam repository I am working on: PM me your email and I can share.

Philip

[Kellis]

Hello all,

I am back again, unfortunately with more issues. Thanks in advance for any tips or help you may have!

As part of a project I am doing, I need to solve a solid experiencing 3D forces. I have created a simple test case based on the largeStrainCantileverBeam tutorial, but am using the elasticSolidFoam solver because I don't expect large displacements in the project. However, moving to 3-D has caused some stability issues which I can't seem to resolve.

Starting from the tutorial case, several changes were made. The applied force was reduced, but acts in all three directions now (the force vector is (1000 1000 1000) on the right patch). Additionally, I changed the blockMeshDict and boundary files so that there are five cells in the third dimension, and the frontAndBack patch is no longer empty.

Additionally, I changed the "planeStress" keyword to "no," since the case is no longer 2-D. This is where the problems start. After a few iterations, I get the "complex eigenvalues detected" warning, and eventually a floating point exception. Interestingly, changing the planeStress flag back to "yes" prevents this from happening, but the solution still doesn't reach tolerance in 10,000 iterations.

So, obviously there has been work done in 3D with the Foam Extend solvers before. Can anyone help identify where my issue is? Attached is the test case where I am getting the error. Thank you!

PS - The end goal is to add a centrifugal force and solve in an SRF, so if anyone has done that before please let me know!

[pippo2013]

Hi everybody!


I am working on a FSI problem, where a thin membrane flaps due to sloshing inside a tank. Due to its flexibility the flapping part might have a contact with the tank.
I think I should use FSI in Foam extend 4.0, namely icoFsiElasticNonLinULSolidFoam,
due to large displacement. However, in order to avoid cell excessive distortion or degeneration when the flap approximates to the bottom I should impose a topology change within the dynamic mesh. But I am afraid in order to include a contact model, I also have to compile a new solver.
Can you please give me any hints?



Philip, would be solids4foam useful in this context?


Thanks in advance!!

[erich1016]

Hi all,

I'm working on the thermal stress simulation during cooling process of a solid plate using elasticThermalSolidFoam. But I found an issue with the reference temperature.

When I set the initial temperature of the plate at 800K and started cooling by applying a fixed 300K boundary condition for all the faces, the result of displacement should show contraction. But I found that I always got expansion no matter I tried cooling or heating. The stress, on the other hand, was as expected. I checked the first time step and found that the plate actually got an initial expansion and then started to contract. This expansion should be due to the temperature difference between 800k and 0K. I know I can set T0 in thermalProperties but it did not seem to make any difference after I tried.

Does anyone know how to avoid this problem?

Thanks.

[bigphil]

Hi Duo,

The stress and displacement will depend on the applied mechanical boundary conditions i.e. boundary conditions for U/DU (or D/DD depending on your version). Can you give some details on your geometry and displacement boundary conditions? In addition, what values of the reference stress-free temperature T0 did you try?

Philip

[erich1016]

Hi Phil,

Thanks for your quick reply.

The geometry is just a 6m * 2m * 1m bar. I'm setting all six faces as displacement free by the following setting:

{
type solidTraction;
traction uniform ( 0 0 0 );
pressure uniform 0;
DT k;
value uniform (0 0 0);
}

As for the reference temperature, I tried to set it the same as the initial temperature of the bar, which is 800K. I assume this could give me a stress free initial state.

About the displacement boundary condition, I know that free surface is not a good condition especially when there are different cooling conditions for different faces. I found that I got unreasonable displacement result (as large as 10^7) when I set top and bottom faces fixed at different temperatures. But the stress results are as expected. Do you have any suggestions on how to set a proper displacement boundary condition that will not introduce rigid body motion and stress concentration?

Thanks so much for your help. I'm new to this area. These solid mechanics solvers are really helpful!

[bigphil]

Hi Duo,

Yes it's important to restrain all rigid body motion; to do this you could:

• create a patch with one face and set it to zero displacement: this will remove rigid translation (rigid rotation could still be an issue);
• set the displacement of one internal cell to be zero; once again this will remove rigid translation (rigid rotation could still be an issue); I have included this feature in soldis4foam (PM me if you want to try it out);
• if there are any symmetry planes, make sure to take advantage of them: these remove rigid translation in the direction normal to the symmetry plane.

If rigid rotation is an issue (it may not be) then you could create 2 or 3 small one-face patches and fix their normal displacement with the fixedDisplacementZeroShear boundary condition.

Philip

[erich1016]

Hi Phil,

Yes, rigid body motion can be restrained if proper displacement boundary is applied, but the problem with the reference temperature for stress free state can actually cause some trouble.

I tried one case where all faces of the bar are constrained by zero displacement. The boundary condition for T is fixed low temperature for all faces to simulate cooling. Due to the contraction after cooling, the surface should be in tension. But the result showed compression instead.

So, how to eliminate the initial expansion is still a problem. This also happens when I tried the solidDisplacementFoam in Openfoam. I found that T0 can be changed in elasticThermalSolidFoam, but changing it makes no differences. If T0 can not be used to set the stress free state, what's the point of changing it?

Thanks,
Duo

[bigphil]

Hi Duo,

Yes, of course, the T0 field should have an effect. Can you share your case here?

Philip

[erich1016]

Hi Phil,

Sorry for the late reply. Here are some of the settings of the case.

blockMeshDict

convertToMeters 1;

vertices
(
    (0 0 0)
    (6 0 0)
    (6 2 0)
    (0 2 0)
    (0 0 1)
    (6 0 1)
    (6 2 1)
    (0 2 1)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (120 40 20) simpleGrading (1 1 1)
);

patches
(
    patch left
	(
          (0 3 7 4)
        )
    pacth right
        (
          (1 2 6 5)
        )   
    patch down
    (
        (0 1 5 4)
    )

    patch up
    (
        (2 3 7 6)
    )
    patch front
    (
        (0 1 2 3)
    )

    patch back
    (
        (4 5 6 7)
    )

);

T

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 800;

boundaryField
{
    left
    {
        type            fixedValue;
        value           uniform 100;
    }
    right
    {
        type            fixedValue;
        value           uniform 100;
    }

    down
    {
        type            fixedValue;
        value           uniform 100;
    }
    up
    {
        type            fixedValue;
        value           uniform 100;
    }
    front
    {
        type            fixedValue;
        value           uniform 100;
    }
    back
    {
        type            fixedValue;
        value           uniform 100;
    }
}

U

dimensions      [0 1 0 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    left
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        DT              k;
        value           uniform (0 0 0);
    }
    right
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        DT              k;
        value           uniform (0 0 0);
    }
    up
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        DT              k;
        value           uniform (0 0 0);
    }
    down
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        DT              k;
        value           uniform (0 0 0);
    }
    front
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        DT              k;
        value           uniform (0 0 0);
    }
    back
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        DT              k;
        value           uniform (0 0 0);
    }
}

thermalProperties

thermal
{
    type     constant;
    C        C [0 2 -2 -1 0 0 0] 434;
    k        k [1 1 -3 -1 0 0 0] 250;
    alpha    alpha [0 0 0 -1 0 0 0] 2.3e-05;
    T0       T0 [0 0 0 1 0 0 0] 800;
}

You can see that I already changed T0 in thermalProperties to 800K.
After running, the displacement results showed expansion in all x,y,z direction. (Sorry, I haven't figure out how to upload images)

Thanks again for your help.

Duo

[bigphil]

Hi Duo,

Can you attach your case as a tgz/zip?

Thanks,
Philip