deleting some liquid bridges in Law2_ScGeom_CapillaryPhys_Capillarity

Capillary pressure is unique in the current code (Law2_ScGeom_CapillaryPhys_Capillarity::CapillaryPressure).
In order to make it different for different contacts, it would need to introduce per-bridge capillarity (CapillaryPhys::CapillaryPressure) which would replace the value from the law in the equations.
It is not very difficult to implement, you could try it.

I was thinking about a "broken" flag, that automatically sets F_cap = 0.
This new flag can be set for every contact/interaction.
I looked in the code, but I need a jump start ...

line 145 in Law2 ... Capillarity.cpp:

Vector3r Fcap = (0.0,0.0,0.0); /// dont if this is correct syntax, but you know what i mean ;)
if (!broken) {
  Fcap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
  }

Oh ok, you don't need different succion in different places, you just want to delete some of the capillary bridges.
In that case, yes, a "broken" flag will do the trick. A good place to use it is line 133 (so that deleting broken interactions and other things will done for you gracefully):

Real Pinterpol = isBroken ? 0 : CapillaryPressure*(R2/liquidTension);

Ok, I solved it in line 145:

Vector3r Fcap;
Fcap = {0.0,0.0,0.0};
 if (!broken) {
 Fcap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
 }

compilation was successfull, but i did not check until now.
Is this ok? Or is your proposal better?

Note sure...

Your change will have a few unpleasant side effects because MeniscusParameters will still be computed with the global suction. For instance, the total fluid volume will include the volume of broken bridges.
Your broken bridges will still have ContactPhysics->meniscus = true, which is not very consistent.

By the way, what is supposed to happen if a new contact is created in the zone where you deleted bridges? For the moment it will create a new bridge since fresh interactions will always have isBroken=false.

Ok, so I will try your way...

My simulation is quasi-static, so there are not so much new contacts within a few steps. But yes, your are right, I was also thinking about that problem. I think I will update the broken flag every step (or every 10 steps), like this:

v_iw = 0.001 #velocity of incr. water-level in [m/s]
dt_now = O.dt
time_cube_full = model_height / v_iw
steps_cube_full = int(time_cube_full / dt_now)
t0 = O.time
for c_final in range (1,steps_cube_full):

  O.step()
  time_now = O.time
  water_height = (time_now-t0) * v_iw

  #determine contacts that are below water_height
  #...
  #set broken=True for that contacts
  #...

What is the best way to make a loop over all interactions, i tried this,

 #determine contacts that are below water_height
 for xx in sphID:
  x = O.bodies[xx]
  if ((x) and (isinstance(x.shape,Sphere))):
   for yy in sphID:
    y = O.bodies[yy]
    if ((y) and (isinstance(y.shape,Sphere))):
     if ((xx != yy) and (O.interactions[xx,yy].isReal)):
      i=O.interactions[xx,yy]
      z_int=i.geom.pos[2] #z coordinate of interaction
      if z_int < water_height:
       print z_int

but i get

  File "5-final.py", line 61, in <module>
    if ((xx != yy) and (O.interactions[xx,yy].isReal)):
IndexError: No such interaction

Ok, i found out, that it is much more easy with

 for i in O.interactions:
  z_int=i.geom.contactPoint[2] #z coordinate of interaction
  if ((z_int < water_height) and (i.isReal)):
   print z_int

but i see negative coordinates, how is this possible?
my model is everywhere z > 0 , but i get:

-0.000325541292072
-0.000760969137883
-0.000368990828872
-0.000660870523873
-0.000396151777372
-0.000625752246923
-0.000506016794083
-0.000366152868505
-0.000415623873535
-0.000442005286835
-0.000364966155445

Ok, i found out, that it is much more easy with

 for i in O.interactions:
  z_int=i.geom.contactPoint[2] #z coordinate of interaction
  if ((z_int < water_height) and (i.isReal)):
   print z_int

but i see negative coordinates, how is this possible?
my model is everywhere z > 0 , but i get:

-0.000325541292072
-0.000760969137883
-0.000368990828872
-0.000660870523873
-0.000396151777372
-0.000625752246923
-0.000506016794083
-0.000366152868505
-0.000415623873535
-0.000442005286835
-0.000364966155445

Hi Jakob,

If you want to loop over interactions you can do like this (with Python, right?)

for i in O.interactions:

To check wether the interaction is Real, you can use this

 if i.isReal:

To check wether the interaction is between spheres or not, you can use something like this

    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):

Then, to access bodies properties, you can do this

 y1 = O.bodies[i.id1].state.pos[1];

Hope it can help

Cheers

You are too quick for me...

Good luck

Something is very strange in my model. My periodic cell is from (0,0,0) to (.002,.002.0042). When I look to the model view, every particle is inside the cell (except mirrored particles). When I print z-positions of my spheres, I get very strange values, that do not agree with the model view. Also I get negative z-positions, which should be impossible!

Has anyone an idea, whats wrong here?

On 19/03/12 07:30, Christian Jakob wrote:
> Question #190529 on Yade changed:
> https://answers.launchpad.net/yade/+question/190529
>
> Status: Answered => Open
>
> Christian Jakob is still having a problem:
> Something is very strange in my model. My periodic cell is from (0,0,0)
> to (.002,.002.0042). When I look to the model view, every particle is
> inside the cell (except mirrored particles). When I print z-positions of
> my spheres, I get very strange values, that do not agree with the model
> view. Also I get negative z-positions, which should be impossible!
Hi Christian,
remember that particles position are saved as absolute values in Yade,
not as relative to the main period (let us say). If you want to know
particles position relative to the cell, you need to wrap them into the
cell first.

Chiara
> Has anyone an idea, whats wrong here?
>

Thanks Chiara, so I have to do this for getting right positions of the spheres, right?

x_cl = 0 #lower x coordinate in front
y_cl = 0 #lower y coordinate in front
z_cl = 0 #lower z coordinate in front
x_cu = 0.002 #upper x coordinate in back
y_cu = 0.002 #upper y coordinate in back
z_cu = 0.003 #upper z coordinate in back

x_pos=[]
y_pos=[]
z_pos=[]
for ii in sphID:
 i = O.bodies[ii]
 if ((i) and (isinstance(i.shape,Sphere))):
  pos = O.bodies[ii].state.pos
  if pos[0] < x_cl:
   pos[0] = pos[0] + (x_cu - x_cl)
  if pos[0] > x_cu:
   pos[0] = pos[0] - (x_cu - x_cl)

  if pos[1] < y_cl:
   pos[1] = pos[1] + (y_cu - y_cl)
  if pos[1] > y_cu:
   pos[1] = pos[1] - (y_cu - y_cl)

  if pos[2] < z_cl:
   pos[2] = pos[2] + (z_cu - z_cl)
  if pos[2] > z_cu:
   pos[2] = pos[2] - (z_cu - z_cl)

Is there no easy way to get right positions?

I implemented this isBroken flag in Law2...Capillarity and in the header files of CapillaryPhys.hpp and Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp
But deleting single liquid bridges is still not working, i tried:

for i in O.interactions:
  z_int=i.geom.contactPoint[2] #z coordinate of interaction
  ...
  if ((z_int < water_height) and (i.isReal)):
    i.isBroken=True #delete liquid bridges

if i look into simulation inspection isBroken flag is false. What i can do now is deleting all bridges with

Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.isBroken=True

but how can i delete single bridges from O.interactions ?

UPDATE

I tried

Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.isBroken=True

and its also not working ...

On 20/03/2012 17:30, Christian Jakob wrote:
> Question #190529 on Yade changed:
> https://answers.launchpad.net/yade/+question/190529
>
> Status: Answered => Open
>
> Christian Jakob is still having a problem:
> I implemented this isBroken flag in Law2...Capillarity and in the header files of CapillaryPhys.hpp and Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp
> But deleting single liquid bridges is still not working, i tried:
>
> for i in O.interactions:
> z_int=i.geom.contactPoint[2] #z coordinate of interaction
> ...
> if ((z_int< water_height) and (i.isReal)):
The problem is most likely with the condition you set: (z_int <
water_height)
Does it work if you only said i.isReal?

Chiara
> i.isBroken=True #delete liquid bridges
>
> if i look into simulation inspection isBroken flag is false. What i can
> do now is deleting all bridges with
>
> Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.isBroken=True
>
> but how can i delete single bridges from O.interactions ?
>

Yes it works, I tested it inserting a counter inside the loop. Its not a problem in the script, but in the source code.

>but how can i delete single bridges from O.interactions ?

What is the result of:

  if ((z_int < water_height) and (i.isReal)):
    print i.isBroken
    i.isBroken=True #delete liquid bridges
    print i.isBroken

?

>Is there no easy way to get right positions?

https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Cell.wrapPt

Thank you Bruno, I will try that wrapPt function.
To the other post:

print i.isBroken

AttributeError: 'Interaction' object has no attribute 'isBroken'

I want to provide more information about the source code.

First I included isBroken flag into the header file of Law2...Capillarity, but there I could only delete all liquid bridges at the same time. So I included isBroken in the header files of CapillaryPhys

https://github.com/jakob-ifgt/trunk/blob/master/pkg/dem/CapillaryPhys.hpp

(for the linear model) and in Ip2_FrictMat_FrictMat_MindlinCapillaryPhys

https://github.com/jakob-ifgt/trunk/blob/master/pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp

(for the hertz model). It seems, that this is not right, because now I can not delete a liquid bridge at all.

Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.isBroken=True is also not working...

print Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.isBroken

AttributeError: type object 'Ip2_FrictMat_FrictMat_MindlinCapillaryPhys' has no attribute 'isBroken'

So the question is: In which header should I put isBroken flag?

Thanks jduriez, that solved my question.

i.phys.isBroken is correct syntax, thank you very much!

Ahah, yes of course, i.phys.isBroken...

Bruno, if you have the time, please check my modifications in Law2...Capillarity::checkFusion (line 223 - 301)

https://github.com/jakob-ifgt/trunk/blob/master/pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp

Compilation is ok, but I get a segmentation fault, when running this with fusionDetection=True.

Christian, you should compile yade in this case with debug=1 option.

Then start your script:

yade-bla-bla-bla --debug scriptname.py

After crash you should get some more useful debugging information.

Anton

I found out, that fusionDetection=True is running, when no particle is deleted. If particles are deleted, e.g. using updSph()

def updSph(rad=1.0):
 i = 0
 for ii in sphID:
  tempID = O.bodies.append(utils.sphere([O.bodies[ii].state.pos[0],O.bodies[ii].state.pos[1],O.bodies[ii].state.pos[2]], material=O.bodies[ii].material.id, radius=O.bodies[ii].shape.radius*rad))
  sphID[i] = tempID
  O.bodies.erase(ii)
  i+=1

, then there is a segmentation fault:

me@debian ~/YADE/my-yade-projects/06-converted-PFC-model >/home/me/YADE/YADE-CapTest-Debug/bin/yade-unknown --debug 2-relaxation.py
Welcome to Yade unknown (debug build)
TCP python prompt on localhost:9001, auth cookie `cskuea'
XMLRPC info provider on http://localhost:21001
Running script 2-relaxation.py

radius of particles increased by factor 2.000000, step 1 of 5
SIGSEGV/SIGABRT handler called; gdb batch file is `/tmp/yade-RIsqC1/tmp-0'
GNU gdb (GDB) 7.0.1-debian
Copyright (C) 2009 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law. Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>.
[Thread debugging using libthread_db enabled]
[New Thread 0x7f62865ef700 (LWP 3159)]
[New Thread 0x7f6296f04700 (LWP 3155)]
[New Thread 0x7f6297705700 (LWP 3154)]
0x00007f62bbdf286d in nanosleep () at ../sysdeps/unix/syscall-template.S:82
82 ../sysdeps/unix/syscall-template.S: Datei oder Verzeichnis nicht gefunden.
        in ../sysdeps/unix/syscall-template.S
Current language: auto
The current source language is "auto; currently asm".
Thread ID 0 not known.

Thread 4 (Thread 0x7f6297705700 (LWP 3154)):
#0 0x00007f62bacb51a3 in select () at ../sysdeps/unix/syscall-template.S:82
#1 0x00000000004f7fd3 in ?? ()
#2 0x00000000004a7ba5 in PyEval_EvalFrameEx ()
#3 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#4 0x00000000004a7752 in PyEval_EvalFrameEx ()
#5 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#6 0x0000000000538a10 in ?? ()
#7 0x000000000041ef47 in PyObject_Call ()
#8 0x0000000000427c1f in ?? ()
#9 0x000000000041ef47 in PyObject_Call ()
#10 0x00000000004a1a53 in PyEval_CallObjectWithKeywords ()
#11 0x00007f629770fcab in ?? () from /usr/lib/pymodules/python2.6/sip.so
#12 0x00007f62979b6498 in ?? () from /usr/lib/pymodules/python2.6/PyQt4/QtCore.so
#13 0x00007f62979e668a in ?? () from /usr/lib/pymodules/python2.6/PyQt4/QtCore.so
#14 0x00007f62b6d60e15 in ?? () from /usr/lib/libQtCore.so.4
#15 0x00007f62bbdea8ca in start_thread (arg=<value optimized out>) at pthread_create.c:300
#16 0x00007f62bacbb92d in clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:112
#17 0x0000000000000000 in ?? ()

Thread 3 (Thread 0x7f6296f04700 (LWP 3155)):
#0 0x00007f62bacb51a3 in select () at ../sysdeps/unix/syscall-template.S:82
#1 0x00000000004f7fd3 in ?? ()
#2 0x00000000004a7ba5 in PyEval_EvalFrameEx ()
#3 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#4 0x00000000004a7752 in PyEval_EvalFrameEx ()
#5 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#6 0x0000000000538a10 in ?? ()
#7 0x000000000041ef47 in PyObject_Call ()
#8 0x0000000000427c1f in ?? ()
#9 0x000000000041ef47 in PyObject_Call ()
#10 0x00000000004a1a53 in PyEval_CallObjectWithKeywords ()
#11 0x00007f629770fcab in ?? () from /usr/lib/pymodules/python2.6/sip.so
#12 0x00007f62979b6498 in ?? () from /usr/lib/pymodules/python2.6/PyQt4/QtCore.so
#13 0x00007f62979e668a in ?? () from /usr/lib/pymodules/python2.6/PyQt4/QtCore.so
#14 0x00007f62b6d60e15 in ?? () from /usr/lib/libQtCore.so.4
#15 0x00007f62bbdea8ca in start_thread (arg=<value optimized out>) at pthread_create.c:300
#16 0x00007f62bacbb92d in clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:112
#17 0x0000000000000000 in ?? ()

Thread 2 (Thread 0x7f62865ef700 (LWP 3159)):
#0 0x00007f62bac8c63d in __libc_waitpid (pid=3160, stat_loc=<value optimized out>, options=0) at ../sysdeps/unix/sysv/linux/waitpid.c:41
#1 0x00007f62bac2a2c9 in do_system (line=<value optimized out>) at ../sysdeps/posix/system.c:149
#2 0x00007f62bac2a600 in __libc_system (line=<value optimized out>) at ../sysdeps/posix/system.c:190
#3 0x00007f62b9eff88a in crashHandler (sig=11) at core/main/pyboot.cpp:45
#4 <signal handler called>
#5 0x00007f62a0fd807b in Law2_ScGeom_CapillaryPhys_Capillarity::checkFusion (this=0x2677240) at /home/me/YADE/trunk/pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp:295
#6 0x00007f62a0fd74b7 in Law2_ScGeom_CapillaryPhys_Capillarity::action (this=0x2677240) at /home/me/YADE/trunk/pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp:178
#7 0x00007f62b16bb3a5 in Scene::moveToNextTimeStep (this=0x26725b0) at /home/me/YADE/trunk/core/Scene.cpp:97
#8 0x00007f62b16c28e5 in SimulationFlow::singleAction (this=0xe46550) at /home/me/YADE/trunk/core/SimulationFlow.cpp:20
#9 0x00007f62b16c3e7e in ThreadWorker::callSingleAction (this=0xe46550) at /home/me/YADE/trunk/core/ThreadWorker.cpp:71
#10 0x00007f62b16c33e3 in ThreadRunner::call (this=0x26728b0) at /home/me/YADE/trunk/core/ThreadRunner.cpp:53
#11 0x00007f62b16c3125 in ThreadRunner::run (this=0x26728b0) at /home/me/YADE/trunk/core/ThreadRunner.cpp:27
#12 0x00007f62b179823d in boost::_mfi::mf0<void, ThreadRunner>::operator() (this=0x2f596e8, p=0x26728b0) at /usr/include/boost/bind/mem_fn_template.hpp:49
#13 0x00007f62b1787f18 in boost::_bi::list1<boost::_bi::value<ThreadRunner*> >::operator()<boost::_mfi::mf0<void, ThreadRunner>, boost::_bi::list0> (this=0x2f596f8, f=..., a=...) at /usr/include/boost/bind/bind.hpp:253
#14 0x00007f62b1776f3b in boost::_bi::bind_t<void, boost::_mfi::mf0<void, ThreadRunner>, boost::_bi::list1<boost::_bi::value<ThreadRunner*> > >::operator() (this=0x2f596e8) at /usr/include/boost/bind/bind_template.hpp:20
#15 0x00007f62b1763284 in boost::detail::function::void_function_obj_invoker0<boost::_bi::bind_t<void, boost::_mfi::mf0<void, ThreadRunner>, boost::_bi::list1<boost::_bi::value<ThreadRunner*> > >, void>::invoke (function_obj_ptr=...) at /usr/include/boost/function/function_template.hpp:153
#16 0x00007f62b183eef3 in boost::function0<void>::operator() (this=0x2f596e0) at /usr/include/boost/function/function_template.hpp:1013
#17 0x00007f62b18179f2 in boost::detail::thread_data<boost::function0<void> >::run (this=0x2f595b0) at /usr/include/boost/thread/detail/thread.hpp:56
#18 0x00007f62b9621200 in thread_proxy () from /usr/lib/libboost_thread.so.1.42.0
#19 0x00007f62bbdea8ca in start_thread (arg=<value optimized out>) at pthread_create.c:300
#20 0x00007f62bacbb92d in clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:112
#21 0x0000000000000000 in ?? ()

Thread 1 (Thread 0x7f62bc1f6700 (LWP 3149)):
#0 0x00007f62bbdf286d in nanosleep () at ../sysdeps/unix/syscall-template.S:82
#1 0x00007f629dce1f66 in pyOmega::wait (this=0xdc6b10) at py/wrapper/yadeWrapper.cpp:370
#2 0x00007f629dce1ac0 in pyOmega::run (this=0xdc6b10, numIter=3000, doWait=true) at py/wrapper/yadeWrapper.cpp:364
#3 0x00007f629ddf8384 in boost::python::detail::invoke<int, void (pyOmega::*)(long, bool), boost::python::arg_from_python<pyOmega&>, boost::python::arg_from_python<long>, boost::python::arg_from_python<bool> > (f=@0xecd828, tc=..., ac0=..., ac1=...) at /usr/include/boost/python/detail/invoke.hpp:94
#4 0x00007f629ddd931a in boost::python::detail::caller_arity<3u>::impl<void (pyOmega::*)(long, bool), boost::python::default_call_policies, boost::mpl::vector4<void, pyOmega&, long, bool> >::operator() (this=0xecd828, args_=0x23f7af0) at /usr/include/boost/python/detail/caller.hpp:223
#5 0x00007f629ddb6549 in boost::python::objects::caller_py_function_impl<boost::python::detail::caller<void (pyOmega::*)(long, bool), boost::python::default_call_policies, boost::mpl::vector4<void, pyOmega&, long, bool> > >::operator() (this=0xecd820, args=0x23f7af0, kw=0x0) at /usr/include/boost/python/object/py_function.hpp:38
#6 0x00007f62b85d9bde in boost::python::objects::function::call(_object*, _object*) const () from /usr/lib/libboost_python-py26.so.1.42.0
#7 0x00007f62b85d9e88 in ?? () from /usr/lib/libboost_python-py26.so.1.42.0
#8 0x00007f62b85e14fb in boost::python::handle_exception_impl(boost::function0<void>) () from /usr/lib/libboost_python-py26.so.1.42.0
#9 0x00007f62b85d6788 in ?? () from /usr/lib/libboost_python-py26.so.1.42.0
#10 0x000000000041ef47 in PyObject_Call ()
#11 0x00000000004a72b8 in PyEval_EvalFrameEx ()
#12 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#13 0x00000000004a9692 in PyEval_EvalCode ()
#14 0x00000000004c98be in PyRun_FileExFlags ()
#15 0x00000000004a0b87 in ?? ()
#16 0x00000000004a7ba5 in PyEval_EvalFrameEx ()
#17 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#18 0x00000000004a7752 in PyEval_EvalFrameEx ()
#19 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#20 0x00000000004a7752 in PyEval_EvalFrameEx ()
#21 0x00000000004a95c1 in PyEval_EvalCodeEx ()
#22 0x00000000004a9692 in PyEval_EvalCode ()
#23 0x00000000004c98be in PyRun_FileExFlags ()
#24 0x00000000004c9ad4 in PyRun_SimpleFileExFlags ()
#25 0x000000000041a6bd in Py_Main ()
#26 0x00007f62bac0ac8d in __libc_start_main (main=<value optimized out>, argc=<value optimized out>, ubp_av=<value optimized out>, init=<value optimized out>, fini=<value optimized out>, rtld_fini=<value optimized out>, stack_end=0x7fffaf4ba2c8) at libc-start.c:228
#27 0x00000000004198d9 in _start ()
Current language: auto
The current source language is "auto; currently c".
Current language: auto
The current source language is "auto; currently asm".
Speicherzugriffsfehler

UPDATE:

segmentation fault also happens, when model is dynamic (when contacts were created and broken rapidly).

Could you open a bug with simple scripts showing the crash?
I'm not really surprised that you can segfault by removing bodies combined with fusion detection.
The 2nd segfault without removing bodies is more surprising, does it happen with older versions (before your last changes)?

> Could you open a bug with simple scripts showing the crash?

I will do so.

> The 2nd segfault without removing bodies is more surprising, does it happen with older versions (before your last changes)?

I cannot test right now, I will make a much more simple script and test this.