Unlogical results in the simulation of Hydrodynamic Interaction for Brownian case

[Nasrollah]

Dear Dr. Ondrej Maxiam,

I am using “SlenderBody” source code to do my simulations. I am modelling the hydrodynamic interaction of sedimentation of a pair of fibres with consideration of Brownian motion. I should mention that I am using periodic boundary conditions for all directions. To do the simulation, I modified four scripts from your source code, including “ThreeShearedFibsBrownian.py”, “RPYVelocityEvaluator.py”, "TemporalIntegrator.py" and “fiberCollection.py”. I have attached them here. Also, I am doing the same simulation, but without considering Brownian motion (ThreeShearedFibsNonBrownian.py), also attached here. The results for non-Brownian case shows that the fibres approach to each other for a period of time. However, for Brownian case, it cannot be seen and they seem to fall down with any interaction on each other. I don’t know if the UAMMD library, which is used for the calculation of hydrodynamic interaction among Brownian particles, works very well or not I would be wondering if you could kindly guide me on this issue.

Best regards, Nasrollah Hajaliakbari Source Code.zip

[omaxian]

Nasrollah,

A few things: 1) I found a small bug in parallelization in the code. I am not ready for you to download my latest updates yet, so for now (since you're only doing 2 filaments), just run with 1 thread (nThr=1). If you want to fix this bug, you can go into FiberCollectionC.cpp and add firstprivate(_MobilityEvaluator) to all of the #pragma omp parallel for statements. 2) The code takes an external force as input, but for gravity you want to pass in a density. You need to add the following methods: a) In FibCollocationDiscretization.py, add

def ForceFromForceDensity(self,ForceDen):
        return np.dot(self._WtildeNx,ForceDen);

b) In fiberCollection.py, add

def ForceFromForceDensity(self,ForceDen):
        RSAns = False;
        if (ForceDen.shape[0] == 3*self._Nfib*self._NXpf):
            ForceDen = np.reshape(ForceDen,(self._Nfib*self._NXpf,3))
            RSAns = True;
        Forces=np.zeros((self._Nfib*self._NXpf,3)); 
        for iFib in range(self._Nfib):
            indsNew = np.arange(iFib*self._NXpf, (iFib+1)*self._NXpf);
            Forces[indsNew,:]=self._fiberDisc.ForceFromForceDensity(ForceDen[indsNew,:]);
        if (RSAns):
            Forces = np.reshape(Forces,3*self._Nfib*self._NXpf);
        return Forces;

c) Then, in the temporal integrator, you need to add in place of line 461,

forceDenExt = self._allFibers.uniformForce([0,0,gravden]);
ForceExt = self._allFibers.ForceFromForceDensity(forceDenExt)

you should check that the sum of the force on each fiber is equal to the density (-10) times the length (1.0). 3) I am not really sure I understand the problem you're describing. If you can send a movie that would be helpful. I can assure you UAMMD is not the problem. Remember, when you include Brownian motion, fibers will rotationally and translationally diffuse, so they won't stay in contact for very long.

Ondrej Maxian

[Nasrollah]

Dear Dr. Maxiam,

Thanks for your quick response and the information about how to properly include density for sedimentation problem. I perhaps was not clear about the hydrodynamic interaction issue – sorry – please let me clarify. We are interested in modelling the inter fibre hydrodynamic interaction (not contact) between two or more fibres falling under gravity, in the case where all fibres are also subject to Brownian dyanamics. Although we do see clear evidence of these interactions without Brownian dynamics, when we turn the Brownian dynamics on, there is no clear evidence of hydrodynamic interactions between fibres, even given the presence of noise. We tried to modify your code to include inter fibre hydrodynamic interactions with Brownian dynamics and have identified a possible issue regarding floating point precision when transferring the data into UAMMD – for instance the comment in RPYVelocityEvaluator.py such as calcBlobTotalVel, calcMOneHalfW, “python will just silently pass by copy and the results will be lost” – we are having a technical problem that we cannot test it. Are we on the right track to get inter fibre hydrodynamic interaction with Brownian motion using your code? Any suggestions you can provide would be gratefully appreciated. If it helps, I have provided two movies for a sedimenting pair of fibres, one with and one without Brownian motion, and whereas inter fibre interactions are evident without Brownian motion, this is not true with Brownian motion. We can provide quantitative results to back this up.

Thank you Nasrollah Hajaliakbari

https://github.com/stochasticHydroTools/SlenderBody/assets/10222516/b8a47a54-7cea-4a33-87d3-3c202538fd48

https://github.com/stochasticHydroTools/SlenderBody/assets/10222516/e2420e43-3de3-4494-9e01-a217f315883c

[omaxian]

Nasrollah,

First of all, you should not have to "modify" any of the source code to get hydrodynamic interactions between the fibers in the Brownian case. All of the code is already there. You should just have to write the main script file. I uploaded an example here which does the falling fibers.

As far as UAMMD, that comment is when you compile UAMMD in the wrong precision. When it's compiled in the wrong precision, you will see it clearly because you will pass in nonzero forces and get back zero velocity. The code prints an error message when that happens. If you're not getting the error message, the velocities should be correct.

When I implemented your example, what I found was that the Brownian simulation agrees with the non-Brownian one for very small values of kbT. As you increase kbT, Brownian motion will overwhelm the hydrodynamic interactions, and you won't "see" the hydrodynamics anymore. My recommendation is that you set kbT to a very small number (like 1e-10) and check that the simulation agrees with the non-Brownian one. Then start increasing kbT and see what happens.

OM

[Nasrollah]

Dear Dr. Ondrej Maxian,

Thanks for your quick response and providing a very useful example. Related to fibre force, I firstly tried to use an approximation. I used the difference between coordinates of each consequent Chebyshev points to approximate the length instead of using the weight matrix to calculate the force on each point. I just wanted to make sure that hydrodynamic interaction works in the source code and then switch that to the accurate option (weight matrix). Anyway, the technical issue was somehow solved and I could run that example you have provided in the source code.

Pertaining to the simulation, I have faced some questions which come as follows. 1- What is this error when I set the parameters in a specific way? “ Traceback (most recent call last): File "/home/home01/scnhv/SlenderBody/Python/Examples/TwoFallingFibs.py", line 143, in maxX,,,_,lamalph,XWithLam = TIntegrator.updateAllFibers(iT,dt,stopcount,Dom,Ewald,gravden=gravden,write=True,outfile=FileString); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/home/home01/scnhv/SlenderBody/Python/TemporalIntegrator.py", line 503, in updateAllFibers lamalph, itsneeded, XWithLam = self.SolveForFiberAlphaLambda(XforNL,XsforNL,iT,dt,tvalSolve,forceExt,lamStar,Dom,Ewald) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/home/home01/scnhv/SlenderBody/Python/TemporalIntegrator.py", line 745, in SolveForFiberAlphaLambda MHalfEta, MMinusHalfEta = self._allFibers.MHalfAndMinusHalfEta(XforNL,Ewald,Dom); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/home/home01/scnhv/SlenderBody/Python/fiberCollection.py", line 1720, in MHalfAndMinusHalfEta MHalfEtaUp = RPYEvaluator.calcMOneHalfW(Xupsampled,Dom); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/home/home01/scnhv/SlenderBody/Python/RPYVelocityEvaluator.py", line 557, in calcMOneHalfW self._GPUEwaldHalf.computeHydrodynamicDisplacements(positions,forces,MHalfW,temperature=0.5,prefactor = 1.0) RuntimeError: [Lanczos] Unknown error (found NaN in result guess) at iteration 10 “.

2- Is there any difference between the use of direct and oversample quadrature for hydrodynamic calculations? How this option must be set?

3- How do NupsampleForDirect is set as a parameter for a simulation?

4- How giters is set for a simulation?

5- Why did you define a new radius? Is this a physics dependent parameter or is it a parameter that relates to the numerical treatment that you have used?

I apologized for the long text. I would be so grateful if you could kindly provide some answers for my questions.

Best regards, Nasrollah Hajaliakbari

[omaxian]

Nasrollah,

0) You should never need to worry about the weights matrix, as all of that is computed by the code internally (at least with the version I updated for you). You only need to specify the force density "grav" as an argument to TIntegrator.updateAllFibers. 1) I have never seen that error before. My guess is that the positions you are using might be going unstable somehow? Does that error happen with smaller time steps? 2) The difference between oversampled and direct quadrature is that the hydrodynamics are computed more accurately with oversampling. When you integrate the RPY kernel, direct quadrature will just directly sum it on the Chebyshev points, while oversampled will actually compute the integrals by oversampling (you can read more about this in Appendix B here). My recommendation is to always use oversampled quadrature.
3) The number of oversampling points should be set based on the fiber aspect ratio. If epsilon=0.004, the number of oversampling points should be approximately 1/epsilon = 250 (you can probably get away with half of this and still get the same results, I would use 100 oversampled points). 4) giters is set for deterministic simulations. For those simulations, there is a time lagging procedure that does the nonlocal velocity in an explicit way without having to solve a big linear system (that is what it does when you set giters=1). Sometimes this process goes unstable and you need to do a few GMRES iterations on the big nonlinear system (this is discussed on pages 27-29 of this pdf). So, in your non-Brownian simulations, set giters=1 until you start to see unstable behavior (I doubt you will for only 2 fibers, but if you do, then increase it). For Brownian simulations, the linear system has to be solved every time step to get the right behavior, so giters has no effect in this case. 5) There are two radii in the problem. There is the physical fiber radius (aspect ratio epsilon), and then there is the corresponding regularized RPY singularity radius epsilon hat. They are related by the formula epsilon hat = exp(3/2)/4*epsilon, which is derived in Appendix A of the paper above. The code will take the real aspect ratio from the user, then convert that to the RPY radius.

Let me know if you have more questions.

OM

[Nasrollah]

Dear Dr. Ondrej Maxian,

Again, thank you for your quick response. At first, based on your previous response before the last one, you mentioned that if I set the kbT to very small values (1e-10), I should expect an agreement between the results of non-Brownian case with Brownian one. I have set the parameters for Brownian case as follows: nFib=2; N=32; Lf=1; nonLocal=1; Ld=5.0; mu=1; eps=1.0/(2109); dt=0.05; omega=0; gam0=0; tf=600.0;
giters = 1; FluctuatingFibs=True; lpstar = 1.0e7; kbT =1.0e-10; Eb=lpstarkbTLf; NupsampleForDirect=200; nThr=16; gravden=-0.1; fibList[0] = DiscretizedFiber(fibDisc,np.reshape(Xs13,3N),np.array([-1.2,0,0])); fibList[1] = DiscretizedFiber(fibDisc,np.reshape(Xs13,3*N),np.array([1.2,0,0]));
For non-Brownian case, I set all parameters the same (Eb=1.0e-3;) However, the results are not same. I can provide a couple of videos to show that if you want.

Next, related to your response for that error, if you set the location of fibres in this way, nFib=2; N=32; Lf=1; nonLocal=1; Ld=50.0; mu=1; eps=1.0/(2109); dt=0.05; omega=0; gam0=0; tf=600.0;
giters = 1; FluctuatingFibs=True; lpstar = 1.0; kbT =1.0e-3; Eb=lpstarkbTLf; NupsampleForDirect=200; nThr=16; gravden=-0.1; fibList[0] = DiscretizedFiber(fibDisc,np.reshape(Xs13,3N),np.array([-1.5,0,0])); fibList[1] = DiscretizedFiber(fibDisc,np.reshape(Xs13,3*N),np.array([1.5,0,0]));

I think you can reproduce the error ((Lanczos] Unknown error (found NaN in result guess)) if you do that. This error can happen even at very small time steps.

Finally, you have created an object called fibDiscFat in fiberCollection.py subroutine. As I went through the source code, I did not really understand why you did this. Can you please provide some explanation for that?

Again, I am sorry for the long text. I would be wondering if you could kindly guide me on these issues.

Kind regards, Nasrollah Hajaliakbari

[omaxian]

Nasrollah,

1) Please ignore the FibDiscFat object for now. I am currently working on a final paper on this code where I will discuss why I am doing that. It is related to using special quadrature for the RPY integrals combined with nonlocal hydro. You are not doing that, since you are doing oversampled quadrature, so this is not relevant for you. I probably should be using a separate branch for this, and I apologize for that. 2) I committed the version of the file I am using to reproduce your results. I ran for 60 seconds instead of 600, and the maximum difference between the Brownian and non-Brownian trajectories was less than 1e-3. That is quite good agreement I think. Obviously when you run longer it will go up (for me it was about 0.01), both because of errors and because you will see translational and rotational thermal noise (real physics). You can get an idea of the time integration error by comparing two different time steps, and an idea of the differences due to noise by comparing two different runs with fluctuations. 3) I was able to reproduce that error in UAMMD. I think the problem is that the Ewald parameter xi (set on line 64, controls the near/far field splitting in the linear scaling algorithm) was too small. This parameter is implementation-dependent and is basically arbitrary. If you increase it, you don't get that error any longer. Note that this is a numerical parameter - so changing it won't affect the results up to numerical error.

Thanks,

Ondrej Maxian

[omaxian]

Nasrollah,

Forgot to mention that you can change two parameters to reduce numerical error: the tolerance for the hydrodynamics (set here) and the tolerance for the iterative solver (set here).

OM

[Nasrollah]

Dear Dr. Ondrej Maxian,

Thank you so much for your reply and providing the info about controlling the residual thresholds in the source code. Also, I am going to analyse the results for this case. To do that, I need to define the Gyration tensor (Ixx, Ixy, Ixz, Iyy, Iyz, Izz) for each fibre. The relation for calculation of this can be found in an article entitled as "Settling dynamics of Brownian chains in viscous fluids" page 3. In that relation, r is the position along the filament, r ̅ is the center of mass of the chain under analysis, and t is time I think the piece of code required for this calculation is this, which must be replaced from line 82 to the end in TwoFallingFibs.py file.

# Time loop
stopcount = int(tf/dt+1e-10);

Gyration=np.zeros((stopcount+1,6*nFib))
X=allFibers.getX()
Xmp=allFibers.getXm()
print('Gyration')
print(X)
print(Xmp)
for iFib in range(nFib):

    for i in range(N): 
        Gyration[0,6*iFib+0]=(X[i+iFib*N,0]-Xmp[iFib,0])*(X[i+iFib*N,0]-Xmp[iFib,0])+Gyration[0,6*iFib+0]  
        Gyration[0,6*iFib+1]=(X[i+iFib*N,0]-Xmp[iFib,0])*(X[i+iFib*N,1]-Xmp[iFib,1])+Gyration[0,6*iFib+1]  
        Gyration[0,6*iFib+2]=(X[i+iFib*N,0]-Xmp[iFib,0])*(X[i+iFib*N,2]-Xmp[iFib,2])+Gyration[0,6*iFib+2]      
        Gyration[0,6*iFib+3]=(X[i+iFib*N,1]-Xmp[iFib,1])*(X[i+iFib*N,1]-Xmp[iFib,1])+Gyration[0,6*iFib+3]    
        Gyration[0,6*iFib+4]=(X[i+iFib*N,1]-Xmp[iFib,1])*(X[i+iFib*N,2]-Xmp[iFib,2])+Gyration[0,6*iFib+4]    
        Gyration[0,6*iFib+5]=(X[i+iFib*N,2]-Xmp[iFib,2])*(X[i+iFib*N,2]-Xmp[iFib,2])+Gyration[0,6*iFib+5]    

    Gyration=(1/Lf)*Gyration    

for iT in range(stopcount): 
    print('Time %f' %(float(iT)*dt));

    maxX,_,_,_ = TIntegrator.updateAllFibers(iT,dt,stopcount,Dom,Ewald,gravden=-0.1,write=True,outfile=FileString);
    X=allFibers.getX()
    Xmp=allFibers.getXm()

    for iFib in range(nFib):

        for i in range(N): 
            Gyration[iT+1,6*iFib+0]=(X[i+iFib*N,0]-Xmp[iFib,0])*(X[i+iFib*N,0]-Xmp[iFib,0])+Gyration[iT+1,6*iFib+0] 
            Gyration[iT+1,6*iFib+1]=(X[i+iFib*N,0]-Xmp[iFib,0])*(X[i+iFib*N,1]-Xmp[iFib,1])+Gyration[iT+1,6*iFib+1]  
            Gyration[iT+1,6*iFib+2]=(X[i+iFib*N,0]-Xmp[iFib,0])*(X[i+iFib*N,2]-Xmp[iFib,2])+Gyration[iT+1,6*iFib+2]      
            Gyration[iT+1,6*iFib+3]=(X[i+iFib*N,1]-Xmp[iFib,1])*(X[i+iFib*N,1]-Xmp[iFib,1])+Gyration[iT+1,6*iFib+3]    
            Gyration[iT+1,6*iFib+4]=(X[i+iFib*N,1]-Xmp[iFib,1])*(X[i+iFib*N,2]-Xmp[iFib,2])+Gyration[iT+1,6*iFib+4]    
            Gyration[iT+1,6*iFib+5]=(X[i+iFib*N,2]-Xmp[iFib,2])*(X[i+iFib*N,2]-Xmp[iFib,2])+Gyration[iT+1,6*iFib+5]   

    Gyration=(1/Lf)*Gyration    
    print(maxX)
np.savetxt('Gyration', Gyration); 

Is it correct or should I modify this? Thanks in advance for any guidance that you would provide.

Kind regards, Nasrollah Hajaliakbari

[omaxian]

Nasrollah,

Here is the correct way to compute that tensor:

# Gyration tensor
X = allFibers.getX();
ChebWts = np.reshape(fibDisc._wX,(fibDisc._Nx,1));
Tensor = np.zeros((3,3,nFib));
for iFib in range(nFib):
    # Compute center of mass
    ThisX = X[iFib*fibDisc._Nx:(iFib+1)*fibDisc._Nx];
    COM = np.dot(ChebWts.T,ThisX)
    Diff = ThisX-COM;
    for iD in range(3):
        for jD in range(3):
            Integrand = Diff[:,iD]*Diff[:,jD]; # element-wise multiplication
            Tensor[iD,jD,iFib]= np.dot(ChebWts.T,Integrand);
print(Tensor[:,:,0]) # This is the first fiber
print(Tensor[:,:,1]) # This is the second fiber

Remember that the midpoint and the center of mass are not the same thing (you need to integrate the positions to get the COM), and that the number of points on each fiber is Nx (not N, the number of tangent vectors). For integration, you need the chebyshev weights w.

OM

[omaxian]

Nasrollah,

One more thing I should mention: because the fibers you are using are very slender (epsilon = 1/2000), you need to use many more oversampling points. You should be using about 1000 points (NupsampleForDirect=1000).

OM

[Nasrollah]

Dear Dr. Ondrej Maxian,

Thank you so much for your reply and guidance.

Best regards, Nasrollah Hajaliakbari