examples/navier_stokes/ex6 compiles, but I get this runtime error when I try and run main2d

examples/navier_stokes/ex6$ ./main2d input2d

+++++++++++++++++++++++++++++++++++++++++++++++++++
At beginning of timestep # 7
Simulation time is 0.7

At end       of timestep # 7
Simulation time is 0.8
+++++++++++++++++++++++++++++++++++++++++++++++++++


Writing visualization files...

P=00000:Program abort called in file ``/home/deleeke/sfw/samrai/2.4.4/linux-g++-debug/include/VisItDataWriter.C'' at line 2945
P=00000:ERROR MESSAGE: 
P=00000:VisItDataWriter<DIM>:
P=00000:    hier::Patch data F_T is less than FLT_MIN 
P=00000:    level: 0  patch: 0  patch_data_id: 46  value: -inf
P=00000:    It cannot be read by VisIt.
P=00000:    Make sure data is properly initialized or
P=00000:    use scale factor to increase its size.
[deleeke-x401a:05263] *** Process received signal ***
[deleeke-x401a:05263] Signal: Aborted (6)
[deleeke-x401a:05263] Signal code:  (-6)

-------------------------------------------------------------

If I want to store
U <- alpha * grad Phi + beta * U

for constant alpha and beta, SC U and CC Phi, one can use HierarchyMathOps::grad() as seen here:
https://github.com/IBAMR/IBAMR/blob/master/src/navier_stokes/StaggeredStokesProjectionPreconditioner.cpp#L322

which works fine because the underlying Fortran code takes into account the grad Phi and beta * U together:
https://github.com/IBAMR/IBAMR/blob/master/ibtk/src/math/fortran/grad2d.f.m4#L430

However, if I want to do the same gradient computation using variable SC alpha, the underlying Fortran code first computes alpha grad Phi terms:
https://github.com/IBAMR/IBAMR/blob/master/ibtk/src/math/fortran/flux2d.f.m4#L328

and then 'axpys' the beta U term:
https://github.com/IBAMR/IBAMR/blob/master/ibtk/src/math/PatchMathOps.cpp#L3585
https://github.com/IBAMR/IBAMR/blob/master/ibtk/src/math/PatchMathOps.cpp#L3644

However in these axpy statements, src2 now holds alphagrad Phi because it is the same patch data as dst, so what is eventually returned by HierarchyMathOps::grad() is
U <- betaalphagrad Phi +alphagrad Phi

instead of

U <- betaU + alphagrad*Phi

==
Similar issues for FC U and src2.
Should we be able to call HierarchyMathOps::grad() with dst == src2 for non-constant alpha (as is the case with constant alpha)?

I'm getting this run-time error in IB/explicit/ex5
This error message is from the main2d, but main3d produces an almost identical error.

This MAY be some issue with my SAMRAI...

/examples/IB/explicit/ex5$ ./main2d input2d

WARNING: Automatically using both VisIt and Silo data writers.
In the future, it will be necessary to include "Silo" in the list of data writer in order to output Silo data.
.
IBStandardInitializer:  Reading from input files.
  base filename: points2d
  assigned to level 0 of the Cartesian grid patch hierarchy


Global seed = 1479860631

IBStandardInitializer:  Deallocating initialization data.

Writing visualization files at timestep # 0 t=0

Writing output files at timestep # 0 t=0
P=00000:Program abort called in file ``/home/deleeke/sfw/samrai/2.4.4/linux-g++-debug/include/Patch.C'' at line 117
P=00000:ERROR MESSAGE: 
P=00000:Failed assertion: d_patch_data[id].isNull()
[deleeke-x401a:01591] *** Process received signal ***
[deleeke-x401a:01591] Signal: Aborted (6)
[deleeke-x401a:01591] Signal code:  (-6)

Ref:

This seems like an unscalable approach. All processors must own atleast one patch on a particular level.
This is not guaranteed in SAMRAI(?).

@knepley Is there a workaround in PETSc?

Currently the Box and Level smoothers for the Stokes + IB equations are virtually identical. The code needs to be refactored. It looks like the Box smoother could be supported by the Level smoother.

In IBFEMethod::computeLagrangianFluidSource(double data_time), the following two lines 888-889 need to move out of the for loop on qp,
// Apply constraints (e.g., enforce periodic boundary conditions) // and add the elemental contributions to the global vector. Q_dof_map.constrain_element_vector(Q_rhs_e, Q_dof_indices); Q_rhs_vec->add_vector(Q_rhs_e, Q_dof_indices);

Ref
https://github.com/IBAMR/IBAMR/blob/f03aed6d0e1fa56d2e03478b1d6346883002db6f/src/IB/IBMethod.cpp#L545-#L560

The matrix d_force_jac should be constructed at each timestep to account for updated position of structure.

if (!d_force_jac) 

should be removed and replaced by

MatDestroy(&d_force_jac);
d_force_jac = NULL;

followed by other existing statements.

Line 173 in LNode-inl.h looks buggy..
it_val = d_node_data_type_arr[k];
It should be
it_val = d_node_data[k];

d_node_data_type_arr contains unique elements and that would always result in a vector of length 1.

This example's main3d compiles and never throws an error, but in never proceeds past this point:

examples/IBFE/explicit/ex4$ ./main3d input3d

 ./main3d input3d*** Warning, This code is deprecated, and likely to be removed in future library versions! ./include/libmesh/auto_ptr.h, line 270, compiled Aug  2 2016 at 11:58:24 ***

IBFEMethod: using FIRST order LAGRANGE finite elements.

*** Warning, This code is deprecated, and likely to be removed in future library versions! ../LIBMESH/src/mesh/boundary_info.C, line 1016, compiled Aug  2 2016 at 12:04:42 ***

To be able to do many-many bodies with CIBMethod this needs to be done programmatically instead of reading vertex files. Folks need this feature to study hydrodynamic suspensions.

Hi, I am doing a project and wondering how to get the distance field if we apply FSM on a mesh that has diffrent refinement level. I saw there are methods about FSM but it's empty.

@boyceg

I am going to enable nightly builds of master on Jenkins. There two options about how the build status is displayed on GitHub:

1. Set build status with a checkmark, red "x", or yellow "unstable/pending" ball on the branches page next to branch name

2. Put build status badge in README

And/or there is the option to get emails:

3. Receive an email only when it is unstable with the console output in the email and another notification then when it returns to passing.

4. Receive an email for every build.

Of those four options which would you like (you can do any combination)?

The present implementation of IBImplicitStaggeredHierarchyIntegrator effectively requires using a StaggeredStokesFACPreconditioner with a StaggeredStokesIBLevelRelaxationFACOperator. The code should be made a bit more general (e.g. using a simple base class for Stokes+IB solvers derived from the StaggeredStokesSolver base class and a new factory class for Stokes+IB solvers).

Ref:

const bool can_be_refined = level_number + 2 < max_levels;

Shouldn't this be:

const bool can_be_refined = level_number + 2 <= max_levels;

P.S: This info is not getting used in the code as of now.

The semi implicit advection diffusion solver does not supply an implementation of applyGradientDetectorSpecialized()

@boyceg: Can .m files be excluded in .clang-format?

If one wanted to run a simulation with IBAMR on a T-shaped domain, but without using IB points to make that non-rectangular domain, would that be easy to do? I’m interested in using the ADR with AMR aspects of IBAMR but I don’t want the time step restrictions associated with the IB wall points.

Issue type

Bug

Affected branches/commits

Branch: master
Commit: fdd3bd1

Description of problem

Looks like there is a compilation error in this example when compiling with GCC

I have not looked into it any further, but this looks like an undefined variable problem.

Compilation error output:

In file included from test_main.cpp:2:0:
example.cpp: In function 'bool run_example(int, char**)':
example.cpp:227:72: error: 'boundary_id' was not declared in this scope
                     if (boundary_info->has_boundary_id(elem,side,1) || boundary_id == boundary_info->has_boundary_id(elem,side,3))
                                                                        ^
In file included from test_main.cpp:2:0:
example.cpp:255:72: error: 'boundary_id' was not declared in this scope
                     if (boundary_info->has_boundary_id(elem,side,1) || boundary_id == boundary_info->has_boundary_id(elem,side,3))
                                                                        ^
Makefile:612: recipe for target 'test2d-test_main.o' failed

Steps to reproduce:

Configure IBAMR to build with gcc based compilers.
Attempt to compile the examples in examples/IBFE/explicit/ex7

I get a runtime error when I try and run the IMP example

/IMP/explicit/ex0$ ./main2d input2d
 
*** Warning, This code is deprecated, and likely to be removed in future library versions! ./include/libmesh/auto_ptr.h, line 270, compiled Aug  2 2016 at 11:58:24 ***
P=00000:Program abort called in file ``../../IBAMR/lib/../src/IB/IMPInitializer.cpp'' at line 234
P=00000:ERROR MESSAGE: 
P=00000:Failed assertion: silo_writer
[deleeke-x401a:03573] *** Process received signal ***
[deleeke-x401a:03573] Signal: Aborted (6)
[deleeke-x401a:03573] Signal code:  (-6)
[deleeke-x401a:03573] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x10330)[0x7f2745432330]
[deleeke-x401a:03573] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37)[0x7f273ff79c37]
[deleeke-x401a:03573] [ 2] /lib/x86_64-linux-gnu/libc.so.6(abort+0x148)[0x7f273ff7d028]
[deleeke-x401a:03573] [ 3] ./main2d[0xb3b6b3]
[deleeke-x401a:03573] [ 4] ./main2d[0xafe514]
[deleeke-x401a:03573] [ 5] ./main2d[0x68eb2f]
[deleeke-x401a:03573] [ 6] ./main2d[0x42d5a3]
[deleeke-x401a:03573] [ 7] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5)[0x7f273ff64f45]
[deleeke-x401a:03573] [ 8] ./main2d[0x42b289]
[deleeke-x401a:03573] *** End of error message ***
Aborted (core dumped)

Let me know if you need any more details!

Users should not have to set ghost cell width stuff in input files (as done in PR #90), but currently there is not a good mechanism for communicating ghost cell width requirements to dynamically constructed linear solvers.

On recent clean build of master, I get a run-time error with the 3d version example (2d runs fine).

$ ./main3d input3d 
*** Warning, This code is deprecated, and likely to be removed in future library versions! ./include/libmesh/auto_ptr.h, line 270, compiled Mar 15 2017 at 12:15:26 ***

IBFEMethod: using FIRST order LAGRANGE finite elements.

*** Warning, This code is deprecated, and likely to be removed in future library versions! ../LIBMESH/src/mesh/boundary_info.C, line 1016, compiled Mar 15 2017 at 12:18:09 ***


Writing visualization files...


+++++++++++++++++++++++++++++++++++++++++++++++++++
At beginning of timestep # 0
Simulation time is 0
Assertion `i < this->n_sides()' failed.
i = 0
this->n_sides() = 0


Stack frames: 14
0: libMesh::print_trace(std::ostream&)
1: libMesh::MacroFunctions::report_error(char const*, int, char const*, char const*)
2: libMesh::Quad4::build_side(unsigned int, bool) const
3: libMesh::Face::build_edge(unsigned int) const
4: ./main3d() [0xa96486]
5: ./main3d() [0xaa303e]
6: ./main3d() [0xaa353d]
7: ./main3d() [0x789d20]
8: ./main3d() [0x7829c3]
9: ./main3d() [0x5dd3bd]
10: ./main3d() [0xa38efe]
11: ./main3d() [0x4c1aa4]
12: __libc_start_main
13: ./main3d() [0x4be939]
[0] ../LIBMESH/src/geom/face_quad4.C, line 130, compiled Mar 15 2017 at 12:17:54
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

@amneetb, does this class work? If not, should we remove it from the library for now?

a) Ref:

  d_vertex_offset[level_number].resize(d_vertex_offset[level_number].size() + 1,
                                         mesh_idx == 0 ? 0 : d_num_vertex[level_number][mesh_idx - 1]);

It looks like d_vertex_offset is just storing the number of vertices of previous mesh registered with this routine, instead of sum of all the number of previously registered mesh vertices.

b) This info is getting used in IMPInitializer::getCanonicalLagrangianIndex() routine and initializeLSiloDataWriter():

I noticed that GNU library was not used in example/ConstraintIB/falling_sphere case. However, we can't make examples/run it without installing GNU library.
Is it possible to fix this issue?

Tks.

With the current version of libMesh C++11 is required and UniquePtr is simply an alias for std::unique_ptr. So calling UniquePtr p(NULL) is not allowed anymore.
For backward compatibility, I believe it could be enough to use the default constructors
UniquePtr p(NULL); -> UniquePtr p;
Any thought?

Ref

This should be

static const int pu_stencil_sz = 2*NDIM;

Otherwise there would be bug in 3d.

Currently, retrieving F_dot to calculate PK1 stress is a non-trivial activity. Passing F_dot to calculate PK1 stress could be useful to model solids that exhibit rate dependent behavior.

It appears that libgslcblas is not being linked to in some cases when checking for libgsl (e.g. using clang on the VM).

Try for example tests/Stokes-IB/ex2 with 4 processors.

The Silo writer should try to respect transient links within a simulation.

The branch boyceg/reduce-memory-use needd the version 3.3 or 3.4 of PETSc and the version 0.9.3
of libMesh while the master branch is already using newer versions. If it is possible, it may be convenient to have the same dependencies in all the branches to avoid the worry of having different versions of the same library.

Ref

The above two lines look suspicious. If d_split_forces is true on line 1309, it won't reach the statements of block L 1351.

Also, what is weak form for external surface forces. Are you including it in the definition of PK1*N or it's treated differently?

Would like to be able to do this for the pressure-like field Phi, computed when stress Normalization is active in IBFE.

@boyceg

So, at the moment, all jobs test every PR when it is made and when new commits are pushed.

Additionally, I have set up some new key phrases for the Jenkins PR jobs to watch out for -- these can be in comments or in commit messages I believe.

retest this please builtin to GHPRB, runs all jobs watching this PR
opt test runs opt job only
skip opt skips opt job, put this in commit message or along with "retest this please"

the rest follow same pattern as above

llvm test
skip llvm
gcc test
skip gcc

Does that sound good?
(note: at the moment, all members of IBAMR organization can trigger these builds)

Support for convective outflow boundary conditions of the form dU/dt + U_INF*dU/dx.

If you try restarting a simulation from an example that does not contain an immersed structure (i.e. navier_stokes/ex3), you will get the following error message:

In the future, it will be necessary to include "Silo" in the list of data writer in order to output Silo data.
.
P=00000:Program abort called in file ``../../../IBAMR/ibtk/lib/../src/lagrangian/LSiloDataWriter.cpp'' at line 2404
P=00000:ERROR MESSAGE:
P=00000:Restart database corresponding to LSiloDataWriter not found in restart file.

The error can be traced back to a call in AppInitializer:
0x000000000081d83a in IBTK::AppInitializer::AppInitializer (this=0x10a1a80, argc=, argv=,
default_log_file_name=...) at ../../../IBAMR/ibtk/lib/../src/utilities/AppInitializer.cpp:263

i.e. this line of code, which looks like it is intended to be deprecated anyway:
https://github.com/IBAMR/IBAMR/blob/master/ibtk/src/utilities/AppInitializer.cpp#L263

==

I understand the purpose of this -- if one doesn't register the Silo writer for an IB case, then this will take care of it. However given that it breaks restarts for purely Eulerian cases, do you think it's safe to remove it at this point in time?

Ref

This should be F += ....

The reason being fluid should experience an opposite reaction. Force on the surface of solid is

f_pressure_solid = -p *n * da
f_pressure_fluid = - f_pressure_solid

f_pressure_fluid should be spread on the Eulerian grid

All the IB/implicit examples compile fine, but the examples fail with the message that they lack any definition of stokes_ib_precond_db required by src/IB/IBImplicitStaggeredHierarchyIntegrator.cpp

The reference I find to that is in line 175 in IBImplicitStaggeredHierarchyIntegrator:

   if (!d_solve_for_position)
    {
        Pointer<Database> stokes_ib_pc_db = input_db->getDatabase("stokes_ib_precond_db");
        d_stokes_solver = new IBImplicitStaggeredStokesSolver(object_name, stokes_ib_pc_db);
        ins_hier_integrator->setStokesSolver(d_stokes_solver);
    }

What are valid entries for stokes_ib_precond_db in the input file?
What would be a sane default in IBImplicitStaggeredHierarchyIntegrator.cpp for this if no preconditioner is provided in the input database?
(The problem is this value is being referenced without any check to see if it exists!)