Error with lubrication law solver

Asked by Léo Sénique on 2019-04-03

Hello,

I'm trying to simulate a viscous flow using a lubrication law, but I'm running into an issue with two of the available resolution methods. I was initially using a modified version of the latest yadedaily release (2018.02b-3905f5add2~bionic), but I have found that this problem still occurs with the original version of this release.
I've simplified the script to try to pinpoint the exact issue:

from yade import pack,plot
import math
import random as rand
import numpy as np
from scipy.stats import gaussian_kde

O.dt=1e-6
slope = math.radians(31)
length = 9
width = 2
height = 2
gravityVector = Vector3(9.81*sin(slope),0.0,-9.81*cos(slope))
O.materials.append(FrictMat(density=2500,frictionAngle=math.radians(35.),poisson=0.2,young=1e7,label='Frict'))
O.materials.append(FrictMat(density=2500,frictionAngle=math.radians(0.),poisson=0.2,young=1e7,label='FrictWall'))

DepBox = O.bodies.append(aabbWalls([(0.,0.,-0.1),(length,width,2*height)],thickness=0.0,material='FrictWall',oversizeFactor=1.0,color=[104/255., 111/255., 140/255.]))
partCloud = pack.SpherePack()
partCloud.makeCloud(minCorner=(length/4.*3.,0.,0.),maxCorner=(length,width,0.8*height), rMean=0.05,num=1000)
spIds = partCloud.toSimulation(material='Frict',color=[205./255,175./255,149./255])

theta_method = 0.55
resolution = 2

O.engines=[
 ForceResetter(),
 InsertionSortCollider(
 [
  Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5),
  Bo1_Box_Aabb(),
  Bo1_Wall_Aabb()
 ],
 verletDist=0.1,
 allowBiggerThanPeriod=False),

 InteractionLoop(
 [
  Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1.5),
  Ig2_Box_Sphere_ScGeom6D()
 ],
 [
  Ip2_FrictMat_FrictMat_LubricationPhys(eta=2,eps=0.001),
 ],
 [
  Law2_ScGeom_ImplicitLubricationPhys(
      activateNormalLubrication=True,
      activateTangencialLubrication=True,
      activateTwistLubrication=True,
      activateRollLubrication=True,
      theta=theta_method,
      resolution=resolution,
      warnedOnce=True,
      maxSubSteps=20,
      SolutionTol=1e-6,
      debug=True)
 ]),

 DragEngine(Rho=1200,Cd=0.47, label='Drag', dead=False),
 NewtonIntegrator(gravity=gravityVector, damping = 0.4, label='newton'),
]
Drag.ids=spIds
Drag.Rho=1200.
Drag.Cd=0.47

This runs fine when using resolution=0 (Iterative exact resolution with substepping) in the lubrication law options.

However, when resolution=1 is used (Newton-Rafson dimentionless resolution), the following error fills the console entirely several times over:
WARN /data/trunk/pkg/dem/Lubrication.cpp:132 NRAdimExp_integrate_u: Max Substepping reach: results may be inconsistant F=nan

A different error appears with resolution=2 (Dichotomy dimentionless resolution), filling the console as well:
WARN /data/trunk/pkg/dem/Lubrication.cpp:199 DichoAdimExp_integrate_u: Max iteration reach: d_left=nan F_left=nan d_right=nan F_right=nan
WARN /data/trunk/pkg/dem/Lubrication.cpp:163 DichoAdimExp_integrate_u: Wrong direction
ERROR /data/trunk/pkg/dem/Lubrication.cpp:175 DichoAdimExp_integrate_u: Initial point problem!! d_left=nan F_left=nan d_right=nan F_right=nan

Is something wrong with my script?
I'm also not entirely sure I understand the difference between those methods and if I should just use the one that works. Is there a difference in precision?

Thanks,

Léo

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-04-04
Last reply:
2019-04-04

Hi,

There are multiple solver for this contact law, but in certain case, I found that some can't give accurate results, just because of the problem's stiffness. The iteratives solvers (dichotomy and Newton-Rafson) are introduced for formulations without exact solutions.

Actually, the last solver (3, dimentionless exact solution with contact prediction) will always give results, without any sub-stepping and/or iterations, and should be used by default. Note that it doesn't implement the theta-method and so, is full-implicit formulation.

Léo Sénique (seniquel) said : #2

Hi William,
If 0 is working, does that mean the system has an exact solution and using 1 and 2 is unnecessary?
How do the methods compare performance-wise?

Hi,

Yes and no. In fact, you have two subsystem: one when the two spheres interacts only with lubrication, and one when you have lubrication and contact.

The "0" method solve the previous system then look if the results is in the right condition. If not, it divide the timestep by two and do it again on the two half timesteps, which sometimes leads to high number of recursions. The "1" and "2" solve directly the whole compound system by iterative methods, and the "3" is like the "0", but it predict the transition, thus, no recursion.

Then , the "3" is most powerful in term of calculation, because it's always a fixed number of computations.

Léo Sénique (seniquel) said : #4

Hi,

Thank you for your answer. Is the "3" method in development? It doesn't seem to exist in the current build of yadedaily and I can't find any reference to it in the documentation. Using resolution=3 gives me the following error:
WARN /home/leo/myYade3/trunk/pkg/dem/Lubrication.cpp:430 go: Nonexistant resolution method. Using exact (0).

Can you help with this problem?

Provide an answer of your own, or ask Léo Sénique for more information if necessary.

To post a message you must log in.