WARNING! rmin>rmax when running TwoPhaseFlowEngine()

Asked by Saeed

Hi There,
When I run the code below, I will be given the Error when arrives to > unsat.initialization()
I guessed the grains out of the box, which I asked in the code to remove that grain, but the grains were not placed outside the box.
Also, using the lower timestep but the problem remains.

Code :

## encoding: utf-8
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
import pylab
from yade import pack
import numpy as np
from numpy import *
import math

## Define Parameters
num_spheres=5000
compFricDegree=30
confiningS=-12500
RhoW=1000
RhoS=2700
youngModulus=75e6
poissonRatio=0.3

## creat a packing with a specific particle size distribution (PSD)
psdSizes=[0.0016,0.002,0.005,0.006,0.0075,0.01,0.015,0.025,0.035,0.05,0.075] #(mm)
psdCumm=[0.0,0.075,0.21,0.23,0.28,0.37,0.48,0.70,0.85,0.89,1.0] #cumulative
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(.0012,0.0012,0.0012) #initial box = 1.2*1.2*1.2 mm
sp.makeCloud(minCorner=mn,maxCorner=mx,num=num_spheres,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=False,seed=True)
sp.psd(bins=150,mass=False)
sp.toSimulation

## create material, which will be used as default
O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,frictionAngle=radians(compFricDegree),density=RhoS,label='spheres'))
O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,frictionAngle=0,density=0,label='frictionless'))

## create walls around the packing
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
 internalCompaction=False, #external (walls) compaction , turn internal compaction off to keep particles sizes constant
 goal1=confiningS,
 goal2=confiningS,
 goal3=confiningS,
 label="triax"
)

newton=NewtonIntegrator(damping=0.4)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 newton
]

O.dt=.1*PWaveTimeStep()
O.dynDt=False

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  e=(triax.porosity)/(1-triax.porosity)
  print ' e: ',e,' Sigma: ',-triax.stress(3)[1]
  if unb<0.15 and abs(confiningS-triax.meanStress)/(-confiningS)<0.015 and e<=1.05:
    break

################
## Start test ###
################
triax.stressMask=2
triax.goal1=triax.goal3=0
goalTop=triax.stress(3)[1]
triax.goal2=goalTop
triax.wall_bottom_activated=False

## Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius)/2.0

## set boundary conditions, the drainage is controlled by decreasing Water pressure and keeping air pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondValue=[0,0,0,0,0,0]
unsat.isPhaseTrapped=False
unsat.defTolerance=0.05
unsat.meshUpdateInterval=200
unsat.useSolver=3
unsat.permeabilityFactor=1
unsat.updateTriangulation=True
unsat.surfaceTension = 0.0728
unsat.initialization()

## Drainage
for pg in arange(1.0e-5,15.0,0.125):
  unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0]
  unsat.invasion()
  unsat.computeCapillaryForce()
  for b in O.bodies:
    O.forces.setPermF(b.id, unsat.fluidForce(b.id))
  while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    Moisture=((RhoW*unsat.getSaturation(False)*triax.porosity*triax.boxVolume)/(RhoS*triax.particlesVolume))*100
    print ' e: ',(triax.porosity)/(1-triax.porosity),' w: ',Moisture,' PC: ',-unsat.bndCondValue[2]
    if unb<0.5:
      break
  if Moisture<18.5:
      break
--------------------------------------------------------------------.
error:
WARNING! rmin>rmax. rmin=2.10136e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.75633e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.02999e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.50292e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.16551e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.67435e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.45035e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.10136e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.27934e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=2.02674e-05 ,rmax=1e-10
WARNING! rmin>rmax. rmin=1.51404e-05 ,rmax=1e-10

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Chao Yuan
Solved:
Last query:
Last reply:
Revision history for this message
Saeed (saeed0668191) said :
#1

Hi there
I think there is a problem with file 'TwoPhaseFlowEngine.cpp' to calculate rmin and rmax :
A portion of this file (engine) that performs radius calculations is listed below:

----------------------
double TwoPhaseFlowEngine::computeMSPRcByPosRadius(const Vector3r& posA, const double& rA, const Vector3r& posB, const double& rB, const Vector3r& posC, const double& rC)
{
    double e[3]; //edges of triangulation
    double g[3]; //gap radius between solid

    e[0] = (posB-posC).norm();
    e[1] = (posC-posA).norm();
    e[2] = (posB-posA).norm();
    g[0] = ((e[0]-rB-rC)>0) ? 0.5*(e[0]-rB-rC):0 ;
    g[1] = ((e[1]-rC-rA)>0) ? 0.5*(e[1]-rC-rA):0 ;
    g[2] = ((e[2]-rA-rB)>0) ? 0.5*(e[2]-rA-rB):0 ;

    double rmin= (std::max(g[0],std::max(g[1],g[2]))==0) ? 1.0e-11:std::max(g[0],std::max(g[1],g[2])) ;
    double rmax= computeEffRcByPosRadius(posA, rA, posB, rB, posC, rC);
    if(rmin>rmax) { cerr<<"WARNING! rmin>rmax. rmin="<<rmin<<" ,rmax="<<rmax<<endl; }

    double deltaForceRMin = computeDeltaForce(posA, rA, posB, rB, posC, rC, rmin);
    double deltaForceRMax = computeDeltaForce(posA, rA, posB, rB, posC, rC, rmax);
    double effPoreRadius;

    if(deltaForceRMin>deltaForceRMax) { effPoreRadius=rmax; }
    else if(deltaForceRMax<0) { effPoreRadius=rmax; }
    else if(deltaForceRMin>0) { effPoreRadius=rmin; }
    else { effPoreRadius=bisection(posA, rA, posB, rB, posC, rC, rmin,rmax); }
    return effPoreRadius;
}
----------------------------------------

Revision history for this message
Robert Caulk (rcaulk) said :
#2

Comparing this script to your other question related to the same error [1], you appear to have decreased your particle size by 2 orders of magnitude. However, you've kept the same relative time step size. I would guess you should decrease the time step as you decrease particle size...

Tough to help you test this because the script you have provided is far from an MWE [2, point 3]. How long does it take to run this script? I quit running it after 15 minutes. Try to widdle it down so you obtain the error with a quick script that we can be run efficiently for trial and error tests. Not to mention, the process of widdling down the script will likely allow you to figure out the problem on your own.

[1]https://answers.launchpad.net/yade/+question/678426
[2]https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Saeed (saeed0668191) said :
#4

hi Robert,
Thank you for the answer you gave me
I lowered the timestep to : O.dt=0.000001*PWaveTimeStep()
and it took a lot of time to run the code but the problem remains and warnings were shown in the terminal.

Can I ignore this Warnings ?

Revision history for this message
Best Chao Yuan (chaoyuan2012) said :
#5

hi,
such warning implies there are some elongated pore throats in the tessellation.
If the grains are not out the box, it might be induced by the local cavities.
I guess 5000 spheres are not enough for your PSD (rmax/rmin>40) to generate a REV sample.

Revision history for this message
Saeed (saeed0668191) said :
#6

Thanks Chao Yuan, that solved my question.