On TunaBrains code, i just noticed, rotating and or translating the paddle in the test does not move the Body, and then the solver fails...

[damian-666]

Just wondering if you implemented it and if it is stable enough to use...

I was thinking to apply drag, from thw wikipedia drag equation , then push the fluid away from the normals of the immersed thing. but i'd rather do it that way CBatty did.. i have some of that some... i can put the link.. its on a big page on Variational methods.. it gives away 2d code for gas ( with 3 circles... i could move those.. and it pushes away the liquid.. with a one frame delay, it seems.. but the shape is nice, and it does not miss any of the fluid... I see code that mentions forward Eulear.. so i was hoping you know something about it.. .. I will report to Tunabrain as well...I think the temperature gradients is the focus of this code.....

[g-amador]

Not sure what you mean!? English appears to have been written to much in a rush...

This code examples use forward explicit Euler. Implicit Euler commonly know as stable fluids from Jos Stam is something different. If you mean viscous drag, viscosity is a parameter so I guess you can alter it to simulate viscous drag. If you mean you want to handle moving solids instead of static ones than the problem that velocities have a different behavior with moving vs static objects: I couldn't find a more simple article (I know I saw once about 7 years ago) but this thesis fairly overviews the boundary conditions: https://www.nada.kth.se/utbildning/grukth/exjobb/rapportlistor/2007/rapporter07/bongart_robert_07018.pdf

[damian-666]

Sorry i have been doing CFD research non stop for 5 months. I am porting your code that your ported from Incremental fluiids, to c#. but before i implement it , I was wondering if the Body ( immersed bodies) can rotate and or move... and push the matter out of the way.. and if you had any issues with it. I could help fix the issues maybe.. I see your version of rotate on a body is different than tunabrains version. thast version does not move the paddle in the example at all, and then the solver fails. I tried a simple angular velocity of 0.1 on the body, in sample 6.. or any sample witha body. ill show you your code to refresh your memory, it is different than the c version...

[damian-666]

Well , it looks like you modified the distance normal.. that could push away stuff. , if I put a angular velocity, i hope that the box pushes the markers and ink around in a vortex... Since it is Euler, I guess i would need to limit how fast i move it. So i think i wlll try it.. unless you tell me it never worked or was too unstable..

@Override public Vector2D closestSurfacePoint(double x, double y) { x -= _posX; y -= _posY; Vector2D result = rotate(x, y, -_theta); x = result.getX(); y = result.getY(); double dx = abs(x) - _scaleX * 0.5; double dy = abs(y) - _scaleY * 0.5; if (dx > dy) { x = nsgn(x) * 0.5 * _scaleX; } else { y = nsgn(y) * 0.5 * _scaleY; } result = rotate(x, y, _theta); x = result.getX() + _posX; y = result.getY() + _posY; return new Vector2D(x, y); } @Override public Vector2D distanceNormal(double nx, double ny, double x, double y) { x -= _posX; y -= _posY; Vector2D result = rotate(x, y, -_theta); x = result.getX(); y = result.getY(); if (abs(x) - _scaleX * 0.5 > abs(y) - _scaleY * 0.5) { nx = nsgn(x); ny = 0.0; } else { nx = 0.0; ny = nsgn(y); } return rotate(nx, ny, _theta); }

On Mon, Sep 12, 2016 at 12:49 AM, Gonçalo Amador notifications@github.com wrote:

Not sure what you mean!? English appears to have been written to much in a rush...

This code examples use forward explicit Euler. Implicit Euler commonly know as stable fluids from Jos Stam is something different. If you mean viscous drag, viscosity is a parameter so I guess you can alter it to simulate viscous drag. If you mean you want to handle moving solids instead of static ones than the problem that velocities have a different behavior with moving vs static objects: I couldn't find a more simple article (I know I saw once about 7 years ago) but this thesis fairly overviews the boundary conditions: https://www.nada.kth.se/utbildning/grukth/exjobb/rapportlistor/2007/ rapporter07/bongart_robert_07018.pdf

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246190282, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmWUJSCmAL9EULtWoqd6-QMjXkMECks5qpDEfgaJpZM4J56kt .

[g-amador]

I just ported the c++ code honestly never attempted to modify it (wanted to port to 3D but haven't got the time since finishing the PhD).

I have prior knowledge with backward Euler in stable fluids from the Masters.

I assume its not just a matter of adding/replacing velocity honestly, it appears it will filter wanted effects.

I would suggest trying to couple two fluids instead, where (the moving solids) density would be tweaked to remain a solid, or if you prefer at each fluid update ignore diffusion of the fluid (solid moving objects).

You could also try to alter locally the velocity in the border cells of the moving object, i.e., the moving obstacle velocity plus the velocity in nearby grid cells, or all together replace the velocity of the cells in the solid border be replaced by the moving obstacle velocity.

[g-amador]

You might also require to add artificial vorticity for this around the moving obstacle.

[damian-666]

i guess by question is... does this code handle moving or rotating boxes. ? it looks different there...more promising..

On Mon, Sep 12, 2016 at 2:36 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

Well , it looks like you modified the distance normal.. that could push away stuff. , if I put a angular velocity, i hope that the box pushes the markers and ink around in a vortex... Since it is Euler, I guess i would need to limit how fast i move it. So i think i wlll try it.. unless you tell me it never worked or was too unstable..

@Override public Vector2D closestSurfacePoint(double x, double y) { x -= _posX; y -= _posY; Vector2D result = rotate(x, y, -_theta); x = result.getX(); y = result.getY(); double dx = abs(x) - _scaleX * 0.5; double dy = abs(y) - _scaleY * 0.5; if (dx > dy) { x = nsgn(x) * 0.5 * _scaleX; } else { y = nsgn(y) * 0.5 * _scaleY; } result = rotate(x, y, _theta); x = result.getX() + _posX; y = result.getY() + _posY; return new Vector2D(x, y); } @Override public Vector2D distanceNormal(double nx, double ny, double x, double y) { x -= _posX; y -= _posY; Vector2D result = rotate(x, y, -_theta); x = result.getX(); y = result.getY(); if (abs(x) - _scaleX * 0.5 > abs(y) - _scaleY * 0.5) { nx = nsgn(x); ny = 0.0; } else { nx = 0.0; ny = nsgn(y); } return rotate(nx, ny, _theta); }

On Mon, Sep 12, 2016 at 12:49 AM, Gonçalo Amador <notifications@github.com

wrote:

Not sure what you mean!? English appears to have been written to much in a rush...

This code examples use forward explicit Euler. Implicit Euler commonly know as stable fluids from Jos Stam is something different. If you mean viscous drag, viscosity is a parameter so I guess you can alter it to simulate viscous drag. If you mean you want to handle moving solids instead of static ones than the problem that velocities have a different behavior with moving vs static objects: I couldn't find a more simple article (I know I saw once about 7 years ago) but this thesis fairly overviews the boundary conditions: https://www.nada.kth.se/utbildning/grukth/exjobb/rapportlist or/2007/rapporter07/bongart_robert_07018.pdf

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246190282, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmWUJSCmAL9EULtWoqd6-QMjXkMECks5qpDEfgaJpZM4J56kt .

[damian-666]

thanks.. that is my plan... i tried it with C Batty variational liquid 2d code ( cant find it now)...andi could push the particles..i will have to interpolate i guess... and see how C Batty code worked.. unless i use flip , then just push the particle...

On Mon, Sep 12, 2016 at 3:08 AM, Gonçalo Amador notifications@github.com wrote:

I just ported the c++ code honestly never attempted to modify it (wanted to port to 3D but haven't got the time since finishing the Ph.D..

I have prior knowledge with backward Euler in stable fluids from the Masters.

I assume its not just a matter of adding velocity honestly. As far as I know I would suggest trying to couple two fluids instead, where (the moving solids) density would be tweaked to remain a solid, or if you prefer at each fluid update ignore diffusion of the fluid (solid moving objects).

You could also try to alter locally the velocity in the border cells of the moving object, i.e., the moving obstacle velocity plus the velocity in nearby grid cells, or all together replace the velocity of the cells in the solid border be replaced by the moving obstacle velocity.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246197951, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmTo346xHNZbjmroakApgOr3E30wsks5qpFHBgaJpZM4J56kt .

[damian-666]

but what trobules me is how the c++ code didnt solve if i put a 0.001 velocity on the box... are you sure you didnt fix it? i woulld not know how to fix this if it happen.. i guess i should test it... i dont even have java working...

On Mon, Sep 12, 2016 at 3:15 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

thanks.. that is my plan... i tried it with C Batty variational liquid 2d code ( cant find it now)...andi could push the particles..i will have to interpolate i guess... and see how C Batty code worked.. unless i use flip , then just push the particle...

On Mon, Sep 12, 2016 at 3:08 AM, Gonçalo Amador notifications@github.com wrote:

I just ported the c++ code honestly never attempted to modify it (wanted to port to 3D but haven't got the time since finishing the Ph.D..

I have prior knowledge with backward Euler in stable fluids from the Masters.

I assume its not just a matter of adding velocity honestly. As far as I know I would suggest trying to couple two fluids instead, where (the moving solids) density would be tweaked to remain a solid, or if you prefer at each fluid update ignore diffusion of the fluid (solid moving objects).

You could also try to alter locally the velocity in the border cells of the moving object, i.e., the moving obstacle velocity plus the velocity in nearby grid cells, or all together replace the velocity of the cells in the solid border be replaced by the moving obstacle velocity.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246197951, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmTo346xHNZbjmroakApgOr3E30wsks5qpFHBgaJpZM4J56kt .

[damian-666]

with the tunabrain code.. i at a angular vel of 0.01 and it blows up .. never converges..

On Mon, Sep 12, 2016 at 3:16 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

but what trobules me is how the c++ code didnt solve if i put a 0.001 velocity on the box... are you sure you didnt fix it? i woulld not know how to fix this if it happen.. i guess i should test it... i dont even have java working...

On Mon, Sep 12, 2016 at 3:15 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

thanks.. that is my plan... i tried it with C Batty variational liquid 2d code ( cant find it now)...andi could push the particles..i will have to interpolate i guess... and see how C Batty code worked.. unless i use flip , then just push the particle...

On Mon, Sep 12, 2016 at 3:08 AM, Gonçalo Amador <notifications@github.com

wrote:

I just ported the c++ code honestly never attempted to modify it (wanted to port to 3D but haven't got the time since finishing the Ph.D..

I have prior knowledge with backward Euler in stable fluids from the Masters.

I assume its not just a matter of adding velocity honestly. As far as I know I would suggest trying to couple two fluids instead, where (the moving solids) density would be tweaked to remain a solid, or if you prefer at each fluid update ignore diffusion of the fluid (solid moving objects).

You could also try to alter locally the velocity in the border cells of the moving object, i.e., the moving obstacle velocity plus the velocity in nearby grid cells, or all together replace the velocity of the cells in the solid border be replaced by the moving obstacle velocity.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246197951, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmTo346xHNZbjmroakApgOr3E30wsks5qpFHBgaJpZM4J56kt .

[g-amador]

For forward Euler you can use this code from the "Fluid Simulation for Computer Graphics"

For forward with obstacles you will need Stable Fluids modified to use obstacles http://www.intpowertechcorp.com/GDC03.pdf

Another old alternative would be using http://www.gamasutra.com/view/feature/1549/practical_fluid_dynamics_part_1.php?print=1

[g-amador]

The only difference I handled in the code was that I had to make sure the right values appeared in some cases because c++ handles out of bound (representable numbers) differently from Java so I slightly altered that in the source code.

I literally tested visually and numerical the C++ vs Java code to address this.

[g-amador]

The C++ code slightly cheats using C++ behavior common in C++ programming but in Java or C# might be a problem with ports sometimes.

[g-amador]

If that code considers rotating or moving obstacles: No!

Moving and rotating obstacles velocity requires you to: 1-Identify the cells in which their borders are within the grid (rough aproximation), e.g.:

---

| || || || S || || || | | || || S || S || S || || | | || S || S || S || S || S || | | || || S || S || S || || | | || || || S || || || | | || || || || || ||__ |

[g-amador]

The object must be within the cells but not all cells are entirely engulfed by the objects in this case accuracy would state use the percentage of the cell occupied (more processing time). Instead use the mean between the cell old velocity and the moving object velocity at the cell (might work).

[damian-666]

thanks... i just went over your thesis.. its in 3d.. very ambitious.. im staying with 2d.. i hope that heat will stir up stuff.. i also have some IB peskin like code... i have seens the stam.. and now i might do the FFT way. after 5 months the FFT way is as good as any multigrid... i use panel an a estimated drag .. and i guess i push the stuff.. .. i may use FLIP scheme after all...i just thikn this one is too slow... but the rest is great... up to mixed density... i hope for stratified dust.... in the air from a volcano... using this solver, with4 types of particles .. ash sand dust... and something.. and with the heat...

On Mon, Sep 12, 2016 at 3:22 AM, Gonçalo Amador notifications@github.com wrote:

For forward Euler you can use this code from the "Fluid Simulation for Computer Graphics"

For forward with obstacles you will need Stable Fluids modified to use obstacles http://www.intpowertechcorp.com/GDC03.pdf

Another old alternative would be using http://www.gamasutra.com/view/feature/1549/practical_fluid_ dynamics_part_1.php?print=1

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246198758, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_Lmc-s_bP8ldjLA7-N4nKwgQWd8ziPks5qpFTegaJpZM4J56kt .

[g-amador]

FFT does not allow moving or either type of obstacles !!!!!

[damian-666]

if i want to run your code.. what is the best way , processing? sorry i never used java outside a browser and now it wont let me..

On Mon, Sep 12, 2016 at 3:31 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

thanks... i just went over your thesis.. its in 3d.. very ambitious.. im staying with 2d.. i hope that heat will stir up stuff.. i also have some IB peskin like code... i have seens the stam.. and now i might do the FFT way. after 5 months the FFT way is as good as any multigrid... i use panel an a estimated drag .. and i guess i push the stuff.. .. i may use FLIP scheme after all...i just thikn this one is too slow... but the rest is great... up to mixed density... i hope for stratified dust.... in the air from a volcano... using this solver, with4 types of particles .. ash sand dust... and something.. and with the heat...

On Mon, Sep 12, 2016 at 3:22 AM, Gonçalo Amador notifications@github.com wrote:

For forward Euler you can use this code from the "Fluid Simulation for Computer Graphics"

For forward with obstacles you will need Stable Fluids modified to use obstacles http://www.intpowertechcorp.com/GDC03.pdf

Another old alternative would be using http://www.gamasutra.com/view/feature/1549/practical_fluid_d ynamics_part_1.php?print=1

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246198758, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_Lmc-s_bP8ldjLA7-N4nKwgQWd8ziPks5qpFTegaJpZM4J56kt .

[g-amador]

Once in the fourier domain there is no turning back :/ That's why so many posterior works over Stam attempt to use Gauss-Seidel Jacobi or Conjugate Gradient.

[g-amador]

The source I provide is a Netbeans project just open an in each file with a main run individually. It will then generate the set of frames you can instead draw the triangles.

[g-amador]

Also cheek if you have not already the "Fluid Simulation for Computer Graphics" its free.

and the incomplete article

http://cg.informatik.uni-freiburg.de/intern/seminar/gridFluids_fluid_flow_for_the_rest_of_us.pdf

[damian-666]

I know what you mean.. its a trick... 1024 x1024 grids .... are too good to be true.. but i found a few papers that say it can be done... as long as the screen is one wave or less... and the borders fake the periodicity. the FFT is used only for the pressure solve. i dont konw how big can the bodies be and i am surprised if its possible then why is AMG still being used... i will have to stuff the papers.. now i am using AMG... and chimera grids... i dont can about accuracy too much..

On Mon, Sep 12, 2016 at 3:34 AM, Gonçalo Amador notifications@github.com wrote:

The source I provide is a Netbeans project just open an in each file with a main run individually. It will then generate the set of frames you can instead draw the triangles.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246199474, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmfKMeJ3WJS3aYhrK7YasGPFiys5Hks5qpFe6gaJpZM4J56kt .

[damian-666]

sorry about my typos... I am not 100% sure FFT can be used , but i have read 3 papers that tell how to do it. and it is for one cycle only... so, 1024 might be doable... wich is more than i get from AMG... are you still doing CFD?

On Mon, Sep 12, 2016 at 3:39 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

I know what you mean.. its a trick... 1024 x1024 grids .... are too good to be true.. but i found a few papers that say it can be done... as long as the screen is one wave or less... and the borders fake the periodicity. the FFT is used only for the pressure solve. i dont konw how big can the bodies be and i am surprised if its possible then why is AMG still being used... i will have to stuff the papers.. now i am using AMG... and chimera grids... i dont can about accuracy too much..

On Mon, Sep 12, 2016 at 3:34 AM, Gonçalo Amador notifications@github.com wrote:

The source I provide is a Netbeans project just open an in each file with a main run individually. It will then generate the set of frames you can instead draw the triangles.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246199474, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmfKMeJ3WJS3aYhrK7YasGPFiys5Hks5qpFe6gaJpZM4J56kt .

[damian-666]

yes.. Peskin himself has mentioned FFT for pressure solv in an immersed boudary paper.. he said any fast poisson solver will work , but he choose to use an FFT based one.. the paper was on... blood flow i have not tried it and he was not specific. i didnt know if the Immersed items needed to be symmetrical .. but i think no... but it must have both sides periodic. and you put boudardies in the off screen area.. or move stuff in the Y direction ... but first i want to see if your box can move the dust... i like tunabrains smoke... with the heat... and i hate c++.. c# has a good math lib now..

On Mon, Sep 12, 2016 at 3:42 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

sorry about my typos... I am not 100% sure FFT can be used , but i have read 3 papers that tell how to do it. and it is for one cycle only... so, 1024 might be doable... wich is more than i get from AMG... are you still doing CFD?

On Mon, Sep 12, 2016 at 3:39 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

I know what you mean.. its a trick... 1024 x1024 grids .... are too good to be true.. but i found a few papers that say it can be done... as long as the screen is one wave or less... and the borders fake the periodicity. the FFT is used only for the pressure solve. i dont konw how big can the bodies be and i am surprised if its possible then why is AMG still being used... i will have to stuff the papers.. now i am using AMG... and chimera grids... i dont can about accuracy too much..

On Mon, Sep 12, 2016 at 3:34 AM, Gonçalo Amador <notifications@github.com

wrote:

The source I provide is a Netbeans project just open an in each file with a main run individually. It will then generate the set of frames you can instead draw the triangles.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246199474, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmfKMeJ3WJS3aYhrK7YasGPFiys5Hks5qpFe6gaJpZM4J56kt .

[damian-666]

thanks... i wil see how you did it :) .. with distance normal.. i have CBattys variational 2d sample.. and i could do a GPU shader based collision all my particles can be checked it they pass through something thin, from the prior pixel... using a line drawn.. and a hittest shader. 1 bit buffer..... wow sorry about my typos, im very tired. this fluid stuff got too much.. im am looking at japanese 6D code.. i went insane...

On Mon, Sep 12, 2016 at 3:31 AM, Gonçalo Amador notifications@github.com wrote:

The object must be within the cells but not all cells are entirely engulfed by the objects in this case accuracy would state use the percentage of the cell occupied (more processing time). Instead use the mean between the cell old velocity and the moving object velocity at the cell (might work).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246199321, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmeoTTa0_ARS-9oaXA9UFp2XC-aI-ks5qpFcTgaJpZM4J56kt .

[g-amador]

Case you didn't figure it out Euler solvers history for CG not for other fields (aeronautics) fluid simulation:

90's there was explicit fluids solvers but could not calculate faster that a 0.001s simulation second, i.e., for a 1 second simulation you needed to do off-line priorly 1000 runs of the fluids solver.

late 90's and initial 2000 Stam said assume the fluid is incompressible, and that it keeps as in particules moving the last direction it went and correct the error to ensure incompressibility and filter numerical error. He first gave FFT (a domain technique) that replaces calculations in R by calculations in the fourier domain, the problem in the fourier domain calculation is way faster but as soon as you are in the FFT you cannot change anything, i.e., no obstacles static or moving. Latter he suggested using successive over relaxation using GS.

Latter some experimented with GPU using shaders or CUDA using Jacabi, Red-black Gauss Seidel, multigrid or Conjugate gradient. These allow obstacles but you have to change velocities per case.

2D will often work for larger reso

[g-amador]

I ended the Masters in 2007/2009 ported code to Java to integrate in a Game Engine I made for the PhD not available yet. I attempted to incorporate either this code as my own and expect to port all to OpenCL or CUDA latter on when I'm working as an hobby.

But not doing any CFD now no time. I have SWE and LBM and particles code working but mostly is 2D and never had the time to report to 3D.

[damian-666]

i love the stam and implicit solve and the shader solvers... im starting to learn the reinmain , weno... and astro code.... everything is connected... but here is the trick .. for FFT you use a twin for each IB you put in.. I am not sure if i can use that trick.. i put it here if you are curious..* it is in bold search this *>>>>>>>>>>>>

Penalty Immersed Boundary Method for an Elastic Boundary with Mass

Yongsam Kim∗ Charles S. Peskin Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012 USA. Preprint submitted to Elsevier Science 26 August 2005 Abstract The immersed boundary (IB) method has been widely applied to problems involv ing a moving elastic boundary that is immersed in fluid and interacting with it. But most of the previous applications of the IB method have involved a massless elastic boundary and used efficient Fourier transform methods for the numerical solutions. Extending the method to cover the case of a massive boundary has required spread ing the boundary mass out onto the fluid grid and then solving the Navier-Stokes equations with a variable mass density. The variable mass density of this previous approach makes Fourier transform methods inapplicable, and requires a multigrid solver. Here we propose a new and simple way to give mass to the elastic boundary and show that the new method can be applied to many problems for which the boundary mass is important. The new method does not spread mass to the fluid grid, retains the use of Fourier transform methodology, and is easy to implement in the context of an existing IB method code for the massless case. Key words: immersed boundary method, penalty method, massive elastic boundary, filament, windsock, flag, inflow velocity, FFT 1991 MSC: 65-04, 65M06, 76D05, 76M20 ∗ Corresponding author. Present address: ICES, University of Texas at Austin, Austin, TX 78712, USA. Email addresses: kimy@ices.utexas.edu (Yongsam Kim), peskin@cims.nyu.edu (Charles S. Peskin). 2 1 Introduction The purpose of this paper is to propose a simple way to extend the immersed boundary (IB) method to a situation in which an elastic moving boundary has mass that cannot be neglected and to show several examples of the application of this new methodology. In the examples chosen, the boundary mass plays a crucial dynamical role. This is in contrast to most previous applications of the IB method, in which the immersed boundary was either massless or neutrally buoyant in the ambient fluid. The IB method was developed to study flow patterns around heart valves, and is a generally useful method for problems in which elastic materials interact with a viscous incompressible fluid. In the IB formulation, the influence of the elastic material immersed in the fluid appears as a localized body force acting on the fluid. This body force arises from the elastic stresses of the material. Moreover, the immersed material is required to move at the local fluid velocity as a consequence of the no-slip condition. The central idea of the IB method is that the Navier-Stokes solver does not need to know anything about the complicated time-dependent geometry of the elastic boundary, and that therefore we can escape from the difficulties caused by the interaction between the elastic boundary and the fluid flow. This whole approach has been applied successfully to problems of blood flow in the heart [18,20–23,25], wave propagation in the cochlea [3,6], platelet aggregation during blood clotting [5], and several other problems [2,8,9,11,13]. In these applications, however, the elastic boundaries have no mass. (In the case of a thick “boundary” like the muscular heart wall, the corresponding 3 assumption is that the boundary is neutrally buoyant, so that the mass of the volume occupied by the boundary is the same as if it were occupied by fluid.) The massless assumption appears in one version of the mathematical derivation of the IB method [19,21,23], see however [7,28], and can be applied when the mass of the moving boundary is too small to have a significant effect on the overall motions of the fluid and the elastic material. With the massless assumption, however, we cannot approach many other problems for which the boundary mass is important. Consider, for example, the simulation of a flapping filament in a flowing soap film as studied by Zhu and Peskin [27,28], who have shown that a massless filament does not flap at all, and that filament mass is essential for the sustained flapping that is seen not only in simulations with mass but also in laboratory experiments [26]. Subsequent experiments (J.Zhang, personal communication) have shown as predicted that an extremely fine filament fails to flap under the same conditions in which flapping of a more massive filament had been seen. Many other examples, in which the massless assumption is not reasonable arise in aerodynamics. Since air is such a light fluid, it is usually the case that elastic boundaries immersed in air have mass that cannot be neglected. Not only the inertial effect of the boundary mass but also the effect of gravity on the boundary mass are important. We study several such examples in this paper. In [7,28], the method used to handle boundary mass was to spread that mass out onto the fluid grid, in much the same manner as boundary force is conven tionally applied to the fluid grid in the IB computations. When this is done, a variable density ρ(x,t) appears in the Navier-Stokes equations, and this com plicates the numerical solver of those equations. Specifically, it makes Fourier methods inapplicable and requires that an iterative method like Multigrid be 4 used instead. This introduces the question of what stopping criterion to use for the iterative solver, i.e., how to optimize the tradeoff between accurate solution of the difference equations at each time step and the computational cost of doing more iterations per time step, a question which does not arise when a direct solution method is available. >>>>>>>>>>>>>>>>>> Our current approach stays much closer to the IB method for the massless case. In particular, it does not spread mass to the fluid and so it requires the same constant-density Navier-Stokes solver as is used in the massless case. The key idea is to introduce a massive boundary point as a twin of each immersed boundary marker where mass is needed. We assume that the massive boundary thus introduced does not interact directly with the fluid and that the original immersed boundary, which will play the same role as in the original IB method, has no mass. The two boundaries are supposed to be the same, and when they move apart, a strong restoring force arises to pull them back together. This is done with a collection of stiff springs connecting the two boundaries. The massive boundary moves according to Newton’s law (F = ma) in which the only forces acting are the forces of the stiff springs and the gravitational force. Meanwhile, at the other end of each spring, the massless boundary marker is moving at the local fluid velocity and spreading force locally to the fluid grid. Among the forces that it spreads is the force on the boundary marker due to the stiff spring (equal and opposite to the force acting on the massive boundary). Because the spring is stiff, the massive boundary point stays close to its twin immersed boundary marker, and the overall effect is that the boundary has mass. Even though the two boundaries (which we call massive and massless bound aries) are supposed to coincide, we allow them to separate by an amount that 5 depends inversely on the stiffness of the spring connecting them. The stiffness parameter of this spring is analogous to a penalty parameter in a constrained optimization problem, which is why we call this new mass-handling idea the penalty immersed boundary (pIB) method. We shall see later that, in the limit of infinite spring stiffness, the two boundaries coincide, and the results obtained by the pIB method approach those of the version of the IB method used by Zhu and Peskin [27,28] for the massive case. 2 Equations of Motion We begin by stating the mathematical formulation of the equations of the motion for a system comprised of a three-dimensional viscous incompressible fluid containing an immersed, elastic boundary with mass [19]. ρ( ∂u ∂t

• u·∇u) =−∇p + µ∆u + f, (1) ∇·u = 0, (2) ρ(x,t) = ρ0 +Z M(r,s)δ(x−X(r,s,t))drds, (3) f(x,t) =Z F(r,s,t)δ(x−X(r,s,t))drds, (4) ∂X ∂t (r,s,t)=u(X(r,s,t),t) =Z u(x,t)δ(x−X(r,s,t))dx, (5) F(r,s,t) = − ∂E ∂X (X(r,s,t),t). (6) Eqs. (1) and (2) are the familiar Navier-Stokes equations for a viscous in compressible fluid subject to an applied body force per unit volume f(x,t). 6 The density ρ(x,t) is the non-uniform mass density of the whole system (fluid+boundary) and µ is the constant fluid viscosity. The unknown func tions in the fluid equations are the fluid velocity, u(x,t); the fluid pressure, p(x,t); and the force per unit volume applied by the immersed boundary to the fluid, f(x,t), where x = (x,y,z) are fixed Cartesian coordinates, and t is the time. Eq. (6) is the immersed boundary equation which is written in Lagrangian form. The parameters r, s label a fixed material point. The unknown X(r,s,t) completely describes the motion of the immersed boundary, and also its spatial configuration at any given time. Note that X(r,s,t) represents a 2-dimensional surface in 3-dimensional space, which is the case in all the 3-D applications in this paper. F = F(r,s,t) is the force density applied by the immersed boundary to the fluid, in the sense that F(r,s,t)drds is the force applied to the fluid by the patch of immersed boundary drds. The elastic contribution to this force density is given by the variational derivative −∂E ∂X of the elastic energy functional E[X(·,·,t)]. This variational derivative is implicitly defined by dE(t) =Z ∂E ∂X (r,s,t)·dX(r,s,t)drds, (7) where dX is a perturbation of the boundary configuration and dE is the resulting perturbation in the elastic energy of the boundary (to first order). In this paper, the energy functional E[X(·,·,t)] is chosen a simple way. In a 2-D space, we consider an immersed boundary as a 1-D rod and use two elasticities: one that resists stretching and compression, and another that resists bending: E[X(·)] = 1 2 csZ (|∂X ∂s |−1)2ds + 1 2 cbZ |∂2X ∂2s | 2 ds, (8) 7 where cs and cb are constants. Regarding a 2-D boundary immersed in a 3 D space as a shell or membrane, one could use quite general strain-stress relations [1,9]. In all our 3-D simulations, however, we use a very simple law to generate the elastic force. The 2-D immersed boundaries are modeled as collections of 1-D rods. Suppose that X(r,s) is an immersed boundary. Then, for a fixed r or s, X(r,s) is a curve which is parameterized as a function of s or r, respectively. Now each such curve can be assigned the same strain-stress relation as is used for the 1-D immersed boundary in our 2-D simulation, and all the forces generated from all these rods are summed up to make an elastic force F(r,s). Eqs. (3), (4) and (5), which we call interaction equations, involve the three dimensional Dirac delta function δ(x) = δ(x)δ(y)δ(z), which expresses the local character of the interaction. ρ0 is the constant fluid density and M(r,s) is the mass density of the immersed boundary in the sense that M(r,s)drds is the mass of the patch drds. In (3) and (4), the delta function δ(x) = δ(x)δ(y)δ(z) is still three-dimensional, but there are only two integrals drds. As a result, ρ(x,t) and f(x,t) are each singular like a one-dimensional delta function with the singularity supported on the immersed boundary. Although ρ(x,t) and f(x,t) are infinite on the immersed boundary, their integrals over any finite volume are finite. In the case of ρ(x,t), for example, we may integrate over an arbitrary region Ω to find Z Ω ρ(x,t)dx=ρ0Vol(Ω) +Z Ω Z M(r,s)δ(x−X(r,s,t))drdsdx =ρ0Vol(Ω) + Z X(r,s)∈Ω M(r,s)drds, (9) 8 which represents the sum of the mass of the immersed material and the fluid that are contained within a fixed space region Ω. Similarly, the integral of f(x,t) over such a volume is the total force applied to the fluid by the part of the immersed material living in the volume. Eq. (5) is the equation of motion of the immersed elastic boundary. It is just the no-slip condition which says that the boundary moves at the local fluid velocity. This is rewritten in terms of the Dirac-delta function in the second form of Eq. (5). We do so in order to expose a certain symmetry with Eq. (4), in which the force generated by the immersed boundary is re-expressed as body force acting on the fluid. This symmetry is important in the construction of our numerical scheme. Note, however, that the integral in Eq. (5) is a triple integral(dxdydz), unlike the integral in Eqs. (3) and (4). Thus ∂X/∂t is finite, unlike ρ and f. Note that Zhu and Peskin [27,28] use this system of equations (1)-(6) and solve it directly. Since the fluid density ρ(x,t) is not a constant, which disables the use of FFT methodology, their numerical solver of the Navier-Stokes equations needs a multigrid method. Now substitute the density defined in Eq. (3) into Eq. (1), and separate the left-hand side of (1) into two terms: one involving the constant density ρ0, and the other containing the singular part of ρ(x,t) that comes from the immersed boundary. Then ρ0( ∂u ∂t
• u·∇u) = −∇p + µ∆u + f −fD, (10) fD(x,t) =Z Du(x,t) Dt M(r,s)δ(x−X(r,s,t))drds, (11) 9 where D Dt is the material derivative. Using the relation: ∂2X ∂t2 (r,s,t) = Du Dt (X(r,s,t),t), (12) we can interpret fD(x,t) in Eq. (11) as an Eulerian body force obtained by using the Dirac delta function to transform the Lagrangian expression FD(r,s,t) = M(r,s) ∂2X ∂t2 (r,s,t), (13) which is known as the D’Alembert force. In a numerical procedure, one might try to evaluate FD(r,s,t) by applying a backward time difference to ∂X/∂t as evaluated by the Eq. (5). Then FD could be applied to the fluid grid like any other immersed boundary force by using the Dirac delta function. This method was proposed by Xiaodong Wang [14] and works when the boundary mass is sufficiently small. Here we use a different way of approximating FD, as described in the following. Now we modify the equations discussed so far to make what we call the penalty immersed boundary (pIB) method [12]. To do this, we introduce two bound aries immersed in a 3-D fluid: one is a massive boundary Y(r,s,t) having all the mass of the elastic immersed boundary, and the other is a massless boundary X(r,s,t). The massive boundary Y(r,s,t) has mass density M(r,s) but does not interact with the surrounding fluid directly. On the other hand, the massless boundary X(r,s,t) interacts with the fluid: it moves with the local fluid velocity and exerts the elastic force locally on the fluid. Both the massless and massive boundaries are supposed to represent the same material surface. If a pair of corresponding boundary points moves apart, a restoring force appears in order to make them move close together (see Fig. 1). The 10 Massive boundary Massless boundary Spring Fig. 1. Massive and massless boundaries are linked together by a collection of a very stiff springs with rest length zero. restoring force FK acting on the massless boundary is FK(r,s,t) = K(Y(r,s,t)−X(r,s,t)), (14) where K is a large constant. Now the massive boundary Y(r,s,t) moves according to the following equa tion: M(r,s) ∂2Y(r,s,t) ∂t2 =−FK(r,s,t) (15) Note that the only force acting on the massive boundary is the reaction force −FK (we shall add to this a gravitational force later), and that the massive boundary does not interact with the fluid directly. On the other end of the spring, the massless boundary X(r,s,t) applies to the nearby fluid both its 11 elastic force FE = −∂E ∂X and the spring force FK. Consider the case in which K in (14) goes to infinity. Then the massive bound ary Y(r,s,t) approaches the massless boundary X(r,s,t), in which case, using X(r,s,t) only, we can see that FK(r,s,t) is same as −FD(r,s,t) of Eq. (13). Note that, in practice, K cannot be infinite but we can keep X(r,s,t) and Y(r,s,t) as close as we like by choosing K sufficiently large. With a large constant K (not infinity), we can use Eqs. (14) and (15) instead of Eq (13) (D’Alembert force) and derive the equations of motion for the pIB method as follows: ρ0( ∂u ∂t
• u·∇u) = −∇p + µ∇2u + f, (16) ∇·u = 0, (17) f(x,t) =Z F(r,s,t)δ(x−X(r,s,t))drds, (18) ∂X ∂t (r,s,t)=u(X(r,s,t),t) =Z u(x,t)δ(x−X(r,s,t))dx, (19) F = FE + FK, (20) FE =− ∂E ∂X , (21) FK(r,s,t) = K(Y(r,s,t)−X(r,s,t)), (22) M(r,s) ∂2Y ∂t2 = −FK(r,s,t)−M(r,s)ge3, (23) where g is the gravitational acceleration and e3 is a unit vector in the vertical direction (against gravity). We have here added a gravitational term on the 12 right-hand side of Eq. (23) in order to allow for the influence of gravity on the massive immersed boundary. (The corresponding term −ρ0ge3 is not included in Eq. (16) because there it only amounts to a redefinition of the pressure.) 3 Numerical Implementation of the pIB method We now describe a formally second-order IB method to solve the equations of motion [15,20]. The word ‘formally’ is used as a reminder that this scheme is only second-order accurate for problems with smooth solutions. Even though our solutions are not smooth (the velocity has jumps in derivative across the immersed boundary), the use of the formally second-order method results in improved accuracy, see [15]. The specific formally second-order method that we use is the one described in [20], generalized here to take into account the massive boundary that is linked to the immersed elastic boundary by stiff springs. In this method, each time step proceeds in two substeps, which are called the preliminary and final substeps. In the preliminary substep, we get data at time level n + 1 2 from data at time level n by a first-order accurate method. Then the final substep starts again at time level n and proceeds to time level n + 1 by a second order accurate method. This Runge-Kutta framework allows the second-order accuracy of the final substep to be the overall accuracy of the scheme. We use a superscript to denote the time level. Thus Xn(r,s) is shorthand for X(r,s,n∆t), where ∆t is the duration of the time step, and similarly for all other variables. Our goal is to compute updated un+1, Vn+1, where V=dY/dt, Xn+1, and Yn+1 from given data un, Vn, Xn and Yn. 13 Before describing how this is done, we have to say a few words about the spatial discretization. There are two such discretizations: one for the fluid and the other for the elastic boundary. The grid on which the fluid variables are defined is a fixed uniform cubic lattice of meshwidth h = ∆x1 = ∆x2 = ∆x3. Now we introduce the central difference operator Di, defined for i = 1,2,3, as follows: (Diφ)(x) = φ(x + hei)−φ(x−hei) 2h , (24) where ei is the unit vector in the i-th coordinate direction. As the notation suggests, the difference operator in i-th direction Di corresponds to the ith component of the differential operator ∇. Thus Dp will be the discrete gradient of p, and D·u will be the discrete divergence of u. We shall also make use the central difference Laplacian L: (Lφ)(x) = 3 X i=1 φ(x + hei) + φ(x−hei)−2φ(x) h2 (25) The fluid mesh and the elastic boundary mesh defined below are connected by a smoothed approximation to the Dirac delta function. It is denoted δh and is of the following form: δh(x) = h−3φ( x1 h )φ( x2 h )φ( x3 h ), (26) 14 where x = (x1,x2,x3), and the function φ is given by φ(r) =                            3−2|r|+√1+4|r|−4r2 8 , if |r|<1 5−2|r|−√−7+12|r|−4r2 8 , if 1≤|r|<2 0 , if 2≤|r| (27) The motivation and derivation for this particular choice is discussed in [21,23]. We are now ready to describe a typical timestep of the numerical scheme. The preliminary substep which goes from time level n to n+ 1 2 proceeds as follows: First, update the position of the massless boundary Xn+ 1 2 (r,s). Xn+ 1 2 −Xn ∆t/2 =X x un(x)δh(x−Xn(r,s))h3 (28) Here and throughout the paper Px denotes the sum over the cubic lattice in physical space on which the fluid variables are defined. Similarly, Pr,s will denote the sum over rectangular lattice in (r,s) space on which the two boundary positions X, Y, the force density F and the mass density M are defined. The key to generalizing the formally second-order method to the massive case is to handle the massive boundary in a manner that closely parallels the numerical treatment of the immersed boundary itself. Thus, we update Yn+ 1 2 in a similar manner to Xn+ 1 2 : Yn+ 1 2 −Yn ∆t/2 = Vn, (29) where Vn is the velocity vector of the massive boundary (the known value at time n∆t, like un). 15 Next, calculate the force density Fn+ 1 2 which is the sum of two parts: one is elastic force and the other is from the spring linked between massless and massive boundaries. F n+ 1 2 K = K(Yn+ 1 2 −Xn+ 1 2 ), (30) Fn+ 1 2 = − ∂E ∂X (Xn+ 1 2 ) + F n+ 1 2 K . (31) These are discretization of Eqs. (20)-(22), respectively. Now we have to change this elastic force density defined on Lagrangian grid points into the force at Eulerian spatial grid points to be applied in the Navier Stokes equations. This is done by a discretization of Eq. (18). fn+ 1 2 (x) =X r,s Fn+ 1 2 (r,s)δh(x−Xn+ 1 2 (r,s))∆s∆r. (32) With fn+ 1 2 in hand, we can turn to solving the Navier-Stokes equations: ρ( u n+ 1 2 i −un i ∆t/2 + 1 2 (u·Dui + D·(uui))n) + Di˜ pn+ 1 2 = µLu n+ 1 2 i + f n+ 1 2 i ,(33) for i = 1,2,3, and D·un+ 1 2 = 0. (34) The notation ˜ pn+ 1 2 is used to distinguish this pressure from the one that is computed when solving Eqs (38)-(39), below. Note that the unknowns in Eqs (33)-(34) are u n+ 1 2 i and ˜ pn+ 1 2 and that they enter into these equations linearly. Since all the coefficients of these equations are constants, the system of Eqs (33)-(34) can be solved by Fast Fourier Transform with the periodic boundary condition [21,23]. 16 Again we have to calculate the velocity Vn+ 1 2 of the massive boundary in the same fashion. M( Vn+ 1 2 −Vn ∆t/2 ) =−F n+ 1 2 K −Mge3. (35) This is nothing but Euler’s method for Eq (23) and completes the preliminary substep. The final substep is the use of un+ 1 2 , Vn+ 1 2 , Xn+ 1 2 and Yn+ 1 2 obtained in the preliminary substep to compute the corresponding quantities at time level n+1 by the second-order accurate midpoint rule. First, using the fluid velocity un+ 1 2 and massive boundary velocity Vn+ 1 2 , we can find the massless boundary configuration Xn+1 and massive boundary position Yn+1. Xn+1 −Xn ∆t =X x un+ 1 2 (x)δh(x−Xn+ 1 2 (r,s))h3, (36) Yn+1 −Yn ∆t = Vn+ 1 2 . (37) The last thing that we have to do is to update the fluid velocity data: ρ( un+1 i −un i ∆t + 1 2 (u·Dui + D·(uui))n+ 1 2 ) + Dipn+ 1 2 = 1 2 µL(un+1 i + un i ) + f n+ 1 2 i ,(38) for i = 1,2,3, and D·un+1 = 0. (39) Just as in the case of Eqs. (33)-(34), Eqs (38)-(39) are a constant-coefficient linear system in the unknowns un+1 i , pn+ 1 2 . This linear system is solved by fast 17 Fourier transform. For the velocity of the massive boundary, we have M( Vn+1 −Vn ∆t ) =−F n+ 1 2 K −Mge3. (40) Since we have now computed un+1, Vn+1, Xn+1 and Yn+1, the timestep is complete. To complete the description of the numerical pIB method, we need to explain the boundary conditions imposed along the edges of our computational do main. In all cases, we use periodic boundary conditions which are very conve nient for the solution of the linear systems of Eqs. (33)-(34) and Eqs. (38)-(39) by fast Fourier transform, but we also break the periodicity in various ways depending on the application. In all of the applications considered here, we need to impose an “inflow veloc ity”. This is done on two (or four) adjacent parallel grid planes (or grid lines in the 2D case). It is very important to use two adjacent planes rather just one, in order to avoid the spatial oscillations that would otherwise propagate throughout the domain via the central-difference structure of our numerical scheme. Although we typically choose either the first or last two grid planes in some coordinate direction, this is of no fundamental significance, since the underlying domain is periodic and all grid points are created equal. The way of driving the “inflow velocity” in this paper is to apply an external force per unit volume equal to f0(x,t) =                α0(u0(t)−u(x,t)) , x ∈ Ω0 0 , otherwise, (41) 18 where Ω0 is the set of grid points comprised of the two grid planes on which we want to control the velocity, u0(t) is the desired velocity on those planes, and α0 is a constant. When α0 is large, the grid velocity is driven rapidly towards u0(t) within Ω0. Note that, when we refine the grid by a factor of 2, since the region Ω0 is fixed, the number of grid planes in Ω0 becomes twice as large. Another way to break the periodicity is to put no-slip walls on faces of the domain. The method that we shall describe for this purpose can indeed be used to put no-slip walls of any shape anywhere within the domain; planar walls along the faces are just a special case. We create fixed no-slip boundaries by laying out an array of “target points” to mark their desired positions. To avoid leaks, the target points should be spaced about half a meshwidth apart (or closer). These target points are fixed in position and do not interact with fluid. Each of them is connected, however, by a stiff spring to a massless immersed boundary point that moves at the local fluid velocity and applies the force generated by the stiff spring locally to the fluid. This provides a feedback mechanism for computing the boundary force needed to enforce the no-slip condition. Note the close analogy between the above construction and the pIB Method, the massive boundary of which is essentially a target position that moves according to Newton’s law of motion instead of being fixed in space. The various methods that we have discussed above for breaking the periodicity of the domain are special cases which hint at the following general remark: In the context of an immersed boundary method, it is not much of a restriction to use periodic boundary conditions for the domain occupied by the fluid. This is because one can immerse within that domain whatever dynamic or static 19 geometry the application demands. As a general setting for immersed bound ary computations, the periodic domain has both philosophical and practical advantages in comparison to other choices one might make. The philosophical advantage is that the periodic domain preserves the translation invariance of free space. Related to this is the practical advantage that we can use the fast Fourier transform to solve the discretized fluid equations at each time step. When walls are wanted in immersed boundary computations, we put them in as immersed boundaries, and this means that a geometrically complicated wall is no harder to model than a planar one. But the boundaries of the fluid domain do not introduce unwanted wall effects, for the simple reason that a periodic domain has, indeed, no boundaries at all. 4 Numerical convergence test of the pIB method In this section, we verify the claim we made in Section 2 that, as K goes to infinity in Eq. (22), the two boundaries X and Y converge to each other and FK=K(X−Y) approaches the negative D’Alembert force −FD. To do that, we consider a unit square of [0 1]×[0 1] filled with an incompressible fluid in which an initially elliptic flexible elastic boundary is immersed. (Throughout this section we use dimensionless variables.) The density and viscosity of the fluid are 1.0 and 0.01, respectively. The initial shape of the elastic boundary is an ellipse with major axis 0.62 and minor axis 0.5. This ellipse is parameterized by arclength x=X0(s), with s in the interval [0,L0], where L0 is the perimeter of the chosen ellipse. (There is no simple formula for an ellipse parameterized by arclength, but the function X0(s) is easily computed numerically.) Of course, s=0 is equivalent to s=L0, 20 since an ellipse is a closed curve. At all times t, the immersed elastic boundary will be described in parametric form x=X(s,t), with X(s,0)=X0(s). Here s is a Lagrangian parameter (a particular value of s corresponds to a particular material point), but s is not equal to arclength in general, since an elastic material by definition may stretch or shrink. It is important to note that the domain of s is fixed: It is always the interval [0,L0] despite any changes in length of the boundary that may occur. This is because the label s attached to a given material point does not change. The immersed boundary with mass density 0.25 interacts with the fluid and generates an elastic force by the following law: F(s) = cs ∂2X ∂s2 (s), (42) where cs=40 is the stiffness coefficient of the immersed boundary. It can be shown that the force law given by Eq. (42) describes a completely flexible boundary (no bending rigidity) in which the tension is equal to cs|∂X/∂s|. This is a special case of Hooke’s law, the special feature being that the rest length of the immersed boundary is zero. Thus, the boundary tries to shrink to a point, but it is prevented from doing so by the incompressibility of the fluid. The equilibrium configuration of the immersed boundary is a circle containing the same area as the initial ellipse. The expected behavior is a damped vibration about this circular state, see Fig. 2. By choosing s=arclength at the initial time, we create uniform tension (equal to cs, since |∂X/∂s|=1) at that time, and in this case the force given by Eq. (42) is normal to the immersed boundary at t=0. This choice of initial condition avoids a rapid transient (boundary layer in time) of tangential equi libration. Since s is not arclength in general, the force given by Eq. (42) may 21 time=0 time=0.064 Fig. 2. Motion of a massive elastic boundary (closed curve). have a nontrivial tangential component at all times other than t=0. Now we choose the penalty spring constant K=200×N2 where N=32, 64, 128, 256, 512, and 1024. For each N, we also change the timestep ∆t=3.2×10−4/N, the space meshwidth ∆x=1/N and the boundary meshwidth ∆s=L0/(4N). That is, when we refine the meshwidths ∆x and ∆s by a factor of 2, the timestep ∆t is also reduced by the same factor. Although this is sufficient for stability within the range of parameters tested, the explicit computation of the boundary force may impose a more severe stability restriction once the parameters have been sufficiently refined, see [24]. Note that the penalty parameter K increases as the meshwidth and timestep are refined. The above manner of increasing it has the property that K∆t2 is constant. In practice, this preserves the numerical stability of the scheme despite the increase in K, as one might expect since the frequency of oscillation of a spring-mass system with spring constant K is proportional to √K. It is fortunate for the overall effectiveness of the pIB method that one can increase K rapidly as the other numerical parameters are refined without inducing numerical instability. The top panel of Fig. 3 shows the maximum distance between massless and 22 0 0.01 0.02 0.03 0.04 0.05 0.06 10 −7 10 −6 10 −5 time maximum distance between massive and massless boundaries N=128 256 512 1024 0 0.01 0.02 0.03 0.04 0.05 0.06 −30 −20 −10 0 10 20 30 40 50 60 time x−directional force acting on a massive boundary point N=128 256 512 1024 Fig. 3. The maximum amplitude ||X − Y||∞ of the oscillation (distance between the two boundaries) and the x-directional force −FK acting on a massive boundary point (‘’ in Fig. 2) are plotted in time for each N. As N gets bigger, the computation gets refined and K increases. As the refinement goes on, the maximum distance of the two boundaries decreases (top panel), and the force acting on the chosen point converges (bottom panel). 23 massive boundaries for each N as a function of time. From the graph (note logarithmic scale, base 10) we can see that the maximum distance of the two boundaries goes to zero as the refinement proceeds, i.e., ||X−Y||∞ → 0 as K →∞, which verifies our first claim. Note in particular that the amplitude of the rapid oscillation in the displacement between the two boundaries also tends to zero as the numerical parameters are refined and K → ∞. The appearance of a constant amplitude of this rapid oscillation is a consequence of the logarithmic scale. The bottom panel of Fig. 3 compares the x-directional force acting on a point of the massive boundary for each N. The chosen point is marked with ‘’ in Fig. 2. The graphs verify our second claim that, K →∞, the force acting on the massive boundary −FK converges (in fact, to FD, but we do not know the exact value FD). Note that the amplitude of the force oscillation appears to converge to zero as the numerical parameters are refined. This is remarkable when one considers that the stiffness parameter is approaching infinity. Overall, Fig. 3 provides empirical evidence that the pIB method works as advertised. One certainly might worry about the high-frequency oscillation introduced into the system by the large value of the penalty stiffness K: In a nonlinear system, such a spurious high-frequency oscillation might interact with the frequency modes of interest in unpredictable ways. But the results in Fig. 3 indicate that this spurious oscillation is not only of high frequency; it is also of low amplitude, in fact, of an amplitude that approaches zero as the numerical parameters are refined. This suggests that the pIB method may indeed be convergent. We provide further empirical evidence for convergence in the remainder of the section. Table 1 shows the convergence ratios of the velocity fields computed by the 24 Table 1 Convergence ratios in L2 and L∞ of the velocity field N ||uN−u2N||2 ||u2N−u4N||2 ||uN−u2N||∞ ||u2N−u4N||∞ ||vN−v2N ||2 ||v2N−v4N ||2 ||vN−v2N||∞ ||v2N−v4N||∞ 32 1.9235 1.1200 2.0367 1.7614 64 2.1044 1.9523 2.0297 1.6366 128 2.0944 1.8466 2.0514 1.8517 256 2.0930 1.8125 2.0430 1.8451 pIB method. Since we do not know the exact solution of the problem, one con vergence ratio needs three numerical solutions for three consecutive N’s. Let (uN,vN) be the velocity field at the chosen time t=0.064, which is roughly half the longest period of the observed vibration of the immersed elastic boundary, and let ||·||2 and ||·||∞ be L2 and L∞ norms, respectively. Table 1 shows that the pIB method has the convergence ratio roughly equal to 2 (at least in the L2 norm; perhaps somewhat worse in the L∞ norm), which implies that the scheme is first order accurate. Although this provides empirical evidence for the convergence of the pIB method, it also underscores the purely formal na ture of its formal second order accuracy. In the formulation of the pIB method (Eqs. (16)-(23)) the appearance of the Dirac delta function induces jumps in certain quantities across the immersed boundary. Specifically, there is a jump in pressure (produced by the normal component of the boundary force) and a jump in the normal derivatives of the velocity components (produced by the tangential component of the boundary force). This lack of smoothness reduces the accuracy of the scheme. These same considerations apply to the original IB method for the massless case just as they do to the pIB method considered here. 25 In this connection, note that Griffith and Peskin [10] have recently described an immersed boundary computation in which the immersed ”boundary” has finite thickness (it is an immersed elastic shell), and in this computation actual second order accuracy is achieved by a formally second-order accurate IB scheme, despite the sudden changes in material properties where the shell meets the fluid. Thus, the case of an immersed elastic structure with finite thickness already has sufficient regularity for the formal accuracy of the IB scheme to manifest itself. Although we have not yet applied the pIB method to immersed elastic structures with finite thickness, it seems reasonable to conjecture that the pIB method of this paper would similarly exhibit actual second-order accuracy in such a case. In the case of a thin immersed boundary, though, it is natural to ask what is the benefit of using a formally second-order accurate IB or pIB scheme, since only first-order accuracy will be realized in practice. We believe that the principal benefit is the reduction in numerical viscosity associated with the skew-symmetric differencing of the nonlinear terms in comparison to a first order upwind scheme. This benefit has in fact been demonstrated in [15], see also [20]. It is not clear, though, whether the skew-symmetric differencing of the nonlinear terms in any sense calls for the second order accurate time dis cretization that is also part of the formal second order accuracy of our scheme. In some applications, it seems possible to use the skew-symmetric differencing of the nonlinear terms together with first-order accurate time discretization, but in other applications this combination appears to be unstable. To show that the reduction of accuracy comes from the spatial discretization and perhaps also from the increase of K, we have also conducted a conver gence study involving time discretization only, with the spatial discretization 26 Table 2 Convergence ratios with the spatial discretization and penalty parameter fixed M ||uM−u2M||2 ||u2M−u4M||2 ||uM−u2M||∞ ||u2M−u4M||∞ ||vM−v2M||2 ||v2M−v4M||2 ||vM−v2M||∞ ||v2M−v4M||∞ 1 4.0075 4.0112 4.0070 4.0070 2 4.0018 4.0028 4.0017 4.0017 4 4.0010 3.9991 4.0009 4.0009 8 3.9934 4.0124 3.9926 3.9926 and the penalty parameter fixed. Table 2 shows the convergence ratios ob tained in this way. We fix N=256 and ∆x=1/N, and change the timestep ∆t=1.25×10−6/M where M=1,2,4,8,16, and 32. The penalty stiffness K=107 is fixed like the spatial discretization. The simulation goes up to time t=0.064. Table 2 shows substantially better convergence ratios than Table 1, and in deed comes to second order accuracy, which would be a convergence ratio of
• In this case, though, convergence is not to the true solution of the partial differential equations, but rather to the solution of the ordinary differential equation system defined by the spatially discretized pIB method, in which time is still continuous, and moreover with constant penalty parameter K. Thus the results shown in Table 2 only verify that our time-stepping scheme is second-order accurate, as claimed. The results in Table 2 have no bearing on the quality of the spatial discretization, nor on the effectiveness of the penalty approach To illustrate the influence of mass on the behavior of an immersed elastic boundary, Figure 4 compares the massive boundary and massless boundary cases by plotting horizontal and vertical radii of the boundary curve in each 27 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 0.62 time radii in horizontal or vertical directions massless case massive case Fig. 4. The influence of boundary mass on the behavior of an immersed elastic boundary in the form of a simple closed curve in a two-dimensional fluid. The curve is initially elliptical, and the horizontal and vertical radii are plotted as functions of time. Both cases show a damped oscillation as the boundary approaches a circular shape, but the massive boundary oscillation has a lower frequency and is sustained longer than that of the massless boundary. The overall downward trend is a form of numerical error (see text) that is somewhat exaggerated here by the vertical scale of the plot. case as functions of time. Of course, the fluid has mass in both cases. In each case, the two radii oscillate in a damped manner, with one radius having its maximum at about the same time as the other radius is at a minimum. As one would expect, the oscillation has a lower frequency and is sustained longer in the massive than in the massless case. The overall downward trend of the plots in Figure 4 represents loss of volume (area), a form of numerical error that is typical of IB computations but can be virtually eliminated by a technique known as ”improved volume conserva 28 tion” [25], which is not implemented here, since most of our applications do not involve enclosed volumes. For other approaches that also lead to better conservation of volume, see [4] and [16]. 5 Applications 5.1 Flapping Filament in a Flowing Soap Film The first application is a flapping flexible filament in a flowing soap film. This replicates a simulation done by Luoding Zhu [27,28] by the method mentioned in the introduction, in which the mass of the immersed boundary is attributed instead to the nearby fluid. Our purpose here is to show that this problem can also be done by the pIB method. The physical experiment that inspired both Zhu’s simulation and ours was performed by Jun Zhang [26] in the Courant Institute WetLab. It involves a vertical soap film that flows downwards between parallel wires under the influence of gravity. A thread, or “filament”, anchored at its upper end but otherwise free, is suspended in the flowing soap film. Two qualitatively dif ferent modes of filament behavior are seen in these experiments. In one, the filament settles into a straight steady configuration, hanging vertically in the flow. In the other, the filament settles into a pattern of sustained flapping, like a flag in the wind. For some choices of parameters, only one of these two behaviors is seen, but there are other choices of parameters for which both behaviors are possible. In the latter (“bistable”) cases, the choice of which behavior occurs is determined by the initial conditions. 29 We have used the same parameters as in [27,28], see Table 3. For the most part, these are also the parameters of the experiment. The most important difference is that the film viscosity has been increased (in both computations) by a factor of 100. The experimental Reynolds number is about 20,000, but the computational Reynolds number is about 200. Before we discuss the results of the filament simulation, we first consider two special aspects of the simulation. The first issue is about how big the spring coefficient K, which appears in Eq. (22), should be to keep two boundaries (massless and massive) close and not to make any computational instability. We choose a constant K by trial and error. As K gets bigger, a larger restoring Table 3 Parameters of flapping filament simulations Parameters magnitude units Film Density 3×10−4 g/cm3 Film Viscosity 1.2×10−3 g/cm·s Film Inflow Velocity (−280)−(−200) cm/s Filament Density 1×10−5 −4×10−4 g/cm Filament rigidity 0.1 erg·cm Filament Length 2.0−3.0 cm Reynolds number 100−210 Gravitational acceleration 980 cm/s2 Domain (rectangle) 9×18 cm×cm 30 force is generated, and the computation becomes more unstable, in which case, we need to reduce the timestep ∆t. Our strategy is to choose an allowed distance between the two boundaries, to adjust K until it is large enough that this allowed distance is not exceeded, and to reduce ∆t, if necessary, as K is increased to avoid numerical instability. Although this process does involve trial and error, the tests for success are extremely simple, and it is easy to zero in on good choices for K and ∆t, as follows: The test for whether K is large enough is simply to monitor the maximum distance between any massless immersed boundary point and its massive twin. Since this distance will be inversely proportional to K (for large K) it is easy to see from the results of a trial run how K needs to be adjusted for the next trial. Numerical instability resulting from too large a choice of ∆t is unmistakable: the immersed boundary seems to explode. Also, in adjusting ∆t when K changes, one can use the rule of thumb that (for a given spatial discretization) stability is determined by K(∆t)2. This has the pleasant consequence that to make K four times larger (for example), one only needs to make ∆t twice as small. Fig. 5 shows that the maximum distance between the two boundaries is controlled within around 1.5 × 10−4 meshwidth. For this simulation, we use K = 108g/cm·s2 and ∆t = 5 × 10−7s. Despite the oscillations seen in the figure, the massless and massive boundaries stay close to each other as the simulation proceeds. The second question is about the way of driving a desired inflow. In section 3 we have suggested a method to drive an inflow which involves an external body force term α0(u0(t)−u(x,t))θ(x) on the right hand side of the Navier-Stokes (momentum) equation. Here, α0 is a large constant, u0(t) = (U0(t),V0(t)) is the desired inflow velocity at time t, and θ(x) is the characteristic function 31 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 −3 4 6 8 10 12 14 x 10 −5 time(s) distance between massive and massless boundaries Fig. 5. The maximum distance in units of meshwidth between the massless and massive boundaries as a function of time. The maximum distance between two boundaries stays below 1.5×10−4 meshwidth here and even after the time shown in this graph. of the chosen region in which we want to specify the velocity. If we use this method, the size of the constant α0 should be taken carefully. In order to make the inflow velocity obtained by solving the Navier-Stokes equations close to the desired inflow in a reasonably short amount of time, the constant α0 should be large. If it is too large, however, α0 can be a source of numerical instability. If such an instability appears, however, it can be eliminated by reducing the time step ∆t. In this way it is possible to control the inflow velocity closely enough for practical purposes. Fig. 6 shows the desired inflow velocity V0(t) and the error of the calculated velocity v(x,t) at one point x in the specified region. The top graph represents the y-component V0(t) of the desired inflow, and the bottom graph is the absolute error of the y-component of the induced velocity v(x,t) from the 32 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 −3 −250 −200 −150 −100 −50 0 time(s) inflow(cm/s) 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x 10 −3 0 1 2 3 4 time(s) abs(V 0 −v(x,t)) Fig. 6. The top panel represents the y-component V0(t) of the desired inflow, and the bottom panel is the absolute error of the y-component of the calculated velocity v(x,t) from the desired velocity V0(t) at one specified point x, i.e., |V0(t)−v(x,t)|. desired velocity. We can see that, except for the very beginning of the time, both velocities match well, i.e. our method of deriving the inflow velocity forces very quickly the velocity in the chosen region to approach the desired inflow. Note that, to avoid a sudden change of the inflow velocity, V0(t) is smoothly increasing in magnitude: v0(1− exp(−t/t0)), where v0 is the final velocity (y-component) and t0 is an arbitrary chosen constant. Also note that, even though this graph shows only a very short time result, the match of both velocities is getting better as the time proceeds. Control of the inflow velocity is accurate to within 1% by the time 3×10−3s have elapsed and continues to improve thereafter. An important result of [27,28] is that filament mass is needed for sustained flapping to occur. Our simulations by the pIB method give the same result. Fig. 7 compares two cases: one has an almost massless filament with density 33 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 time(s) x of free end(cm) mass=0.00001 mass=0.0004 Fig. 7. The x-coordinate of the free ends of two filaments are plotted as functions of time. The sustained flapping motion is observed only in the case of a filament with enough mass. 1×10−5g/cm and the other has a heavier filament with density 4×10−4g/cm, which is the value used in [27,28]. This value is the same as in Zhang’s ex periments [26], allowing for a factor of 2 to account for a bulge of the film that forms around the filament as a result of surface tension. Other conditions are also the same: the inflow velocity is -280cm/s, the filament length is 3cm, and the initial configurations are same as in [27,28]. The figure plots the x coordinate of the free end of the two filaments as functions of time. As time proceeds, while the small mass filament rapidly becomes straight and then re mains straight like a rigid body, the larger mass filament settles into a flapping mode with the total excursion of the free end 1.8cm and the frequency 48/s. Fig. 8 shows the vorticity contours in both cases at two fixed times. Note the concentrated vorticity on the two sides of the domain as well as around the filament and in the wake of the filament. The concentrated vorticity near the side walls represents the boundary layers that are generated there. (Recall 34 Fig. 8. Motion of two filaments with different mass densities: 10−5g/cm (left panels) and 4×10−4g/cm (right panels). The vorticity contours are drawn at two different times: 0.06s (top) and 0.18s (bottom). the manner in which the no-slip condition has been imposed on the side walls through the use of tether points, as explained in Section 3.) The concentrated vorticity on the two sides of the filament similarly represents the boundary layers generated by the shearing effect of the (almost inextensible) filament on the flow. In the case of the more massive filament, discrete vortices of opposite sign are shed as the filament flaps. In order to explore the issue of bistability, we choose a filament of length 2cm and mass density 4×10−4g/cm. The steady inflow velocity v0 is -200cm/s, which produces a Reynolds number of 100. Two different cases are created by two different initial perturbations from the straight configuration of the fila 35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 −0.4 −0.2 0 0.2 0.4 0.6 4 4.5 5 12.5 13 13.5 14 14.5

On Mon, Sep 12, 2016 at 4:05 AM, Gonçalo Amador notifications@github.com wrote:

Case you didn't figure it out Euler solvers history for CG not for other fields (aeronautics) fluid simulation:

90's there was explicit fluids solvers but could not calculate faster that a 0.001s simulation second, i.e., for a 1 second simulation you needed to do off-line priorly 1000 runs of the fluids solver.

late 90's and initial 2000 Stam said assume the fluid is incompressible, and that it keeps as in particules moving the last direction it went and correct the error to ensure incompressibility and filter numerical error. He first gave FFT (a domain technique) that replaces calculations in R by calculations in the fourier domain, the problem in the fourier domain calculation is way faster but as soon as you are in the FFT you cannot change anything, i.e., no obstacles static or moving. Latter he suggested using successive over relaxation using GS.

Latter some experimented with GPU using shaders or CUDA using Jacabi, Red-black Gauss Seidel, multigrid or Conjugate gradient. These allow obstacles but you have to change velocities per case.

2D will often work for larger reso

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246201283, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmSFQh66gqRj2-PDrRtJ4QHJ8Do-Dks5qpF8jgaJpZM4J56kt .

[damian-666]

i am looking at the CIP method.. a japanese guy made a 6D method 20 years ago.. this is for details on the surface flow.. and accurate waves in 2d.. and 3d.. im exhausted though.. adaptive mesh is also the new think.. every week i learn a new scheme and i never get to the bottom.

On Mon, Sep 12, 2016 at 4:08 AM, Gonçalo Amador notifications@github.com wrote:

I ended the Masters in 2007/2009 ported code to Java to integrate in a Game Engine I made for the PhD not available yet. I attempted to incorporate either this code as my own and expect to port all to OpenCL or CUDA latter on when I'm working as an hobby.

But not doing any CFD now no time. I have SWE and LBM and particles code working but mostly is 2D and never had the time to report to 3D.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246201390, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmfPth_FrTzf3aJCm6REN8GEtlHegks5qpF-rgaJpZM4J56kt .

[damian-666]

you might like IVOCK code it will help put the code in the size hat fits in WARP.. it can do cuda and open cl i think..

and people say copy datat to the gpu is costly.. so then they do everything on GPU or use compressed grids and WARP sizes or streams... IVOCK is on github..its a vorton method tho..

On Mon, Sep 12, 2016 at 4:08 AM, Gonçalo Amador notifications@github.com wrote:

I ended the Masters in 2007/2009 ported code to Java to integrate in a Game Engine I made for the PhD not available yet. I attempted to incorporate either this code as my own and expect to port all to OpenCL or CUDA latter on when I'm working as an hobby.

But not doing any CFD now no time. I have SWE and LBM and particles code working but mostly is 2D and never had the time to report to 3D.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246201390, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmfPth_FrTzf3aJCm6REN8GEtlHegks5qpF-rgaJpZM4J56kt .

[g-amador]

Used GPU before with CUDA and OpenCL the times in the thesis actually can be way faster.

[damian-666]

20x faster? or 4x... i read for the fast FFT it doesnt matter.. better do it on the CPU. or where ever the data is ... visual studio update has some texture compression stuff .. and i read if its 64k i think.. its fits in one warp.... i might try Cudaify.. .. i was doing almost realtime mandelbrot at full resolution...with an opencl demo.. so im planning ot use some 1 bit.. 2 bit .. paletted bitmaps..if i have to pass them around.

i would love if openCL finds the processor in skylake, and if the is an AMD or nvidia ,,the we could use them both on newer machines.. skylake GPU on the chip is almost good enough,, not quite ... so alot of systems will have two powerful gpus. i cant seem to query my intel one..but i think its too old..

im still thinkingo of the FFT solv.. there are alot of tricks.. here are a few more papers.. if its just for pressure solv in air...the presence of one small 2d character, would not ruin it..just not be perfectly accurate? the FFT region might just be the sky ... then using some other grids in chimera style....

i guess its pretty easy to try them... for a side scoller game.. the boundaries are often good enough..

http://ms.mcmaster.ca/~bprotas/MATH745c/spectr_01.pdf

http://ufrmeca.univ-lyon1.fr/~buffat/FDMIF/articleIBM.pdf

.https://people.eecs.berkeley.edu/~demmel/cs267_Spr05/Lectur es/Lecture26/Lecture_26_AMR_part2_PhilColella_05.pdf

On Mon, Sep 12, 2016 at 5:06 AM, Gonçalo Amador notifications@github.com wrote:

Used GPU before with CUDA and OpenCL the times in the thesis actually can be way faster.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246204629, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmaSlxuQrAmU6MHvW0vp0Sgat3H8Fks5qpG1RgaJpZM4J56kt .

[damian-666]

Hello again. I have compared your code and yes it appears there is is basically a straight port from tunabrain.

thanks for porting it .. i think i might giveit a port to c# unless you know a better way to do multiple density gas with immersed boundary. I have tried stuff and play to try more with this code and solver. So just letting you know my plan and about a few more things i read about that might interest you.

I have found some but the license is GPL... i dont think a port to c# will clear me. I dont want to show my creature code just because they swim in someone elses fluids..

I think i am going to try it.. i like what i see so far. I do wan air over water eventually.. and definitely Adaptive grids..

but the results for multiple density are good.. when the object is hot.

I dont get convergence on the second pass and so I tried to set the velocity only after the entire solve This is a very weak coupling. Then as you suggested , i move the fluid around the object using a direct setting of the fluid velocity ( will need to be careful of the maximum change or flux allowed) . and apply panel drag to my IB based on the density near the panel.

because my immersed object is small compared to the gas, ( in air) and speed is far more important than accuracy, I tried this with fair results. the gass seems pushed out of the way if the box is thick enough.

It seems, why should the solver need to know the velocity of the IB if it will simply be one frame late? by then it will be reflected in the fluid using a direct forcing as the initial conditions.

-solve pressure. ( without the velociy of IB objects)

-draw

-apply drag.

-apply direct forcing . repeat

I dont know exactly why it didnt converge, i think that the velocity needed to be multipled by DT in some places. However, it still didnt solve.

I find good enough results using the panel method, adding a drag = 1/ vel * vel x panel length normal to the relative flow , at two points on each panel that is a loose body. skipping all small panels.

Ghost matter might be needed to "thicken" panels, or use collision detection for each particle.. if i use shaders to manage the particles and detect collisions using graphics AND code

I dont feel i was to use vorton methods are there are tons..

So, the solver there seems good enough if the grid is less that 256 x256 I will look for a shader version. Or i will try to Cudaify it...

But herer is another IB method i havent explored...

https://github.com/barbagroup/petrbf

I have FLIP code from Martin Quay that uses a graphics Guassian Splat to determine the density graphically. For 2d there are many tricky and "visual" ways to solve for pressure, and some are begging for an all shader solution. for the collision, for the density , the retendering, and the fluid dynamics. I get about 200,000 particles... but its too slow... no temperature.. no level sets...not adaptive grid...

But i am still wondering about using a 1024x1024 grid for the entire atmosphere..

a google for Spectral methods on incompressible fluids shows many ways.

periodic meaning inflow on one end, equal outflow of mass on the other ( and hidden stuff in the margins)

For a two phase free boundary i am not sure. Perhaps two solvers entirely: one for air, using the water surface as an IB.. and taking its velocity in the normal direction at as affecting the fluid velocity near the interface, and vice versa to a different extent based on the mass ratio.

I do see FFT on GPU solutions. https://slashdot.org/story/06/05/29/1424213/high-performance-fft-on-gpus

in python there are many ways to do IB. but i hate python.. i could use Ironpython since it compiles... or .net..

and i find many papers on using th e FFT for a pressure solve.

48x thats what i like. i even see 1000x when using Jacobi free Krylov methods..

http://www.nvidia.com/content/new_cudawmfg/cudabrowser/ assets/data/applications.xml

its really easy to just read and read and read...

just sharing some stuff.

On Mon, Sep 12, 2016 at 9:27 AM, Damian Hallbauer < damian.hallbauer@gmail.com> wrote:

20x faster? or 4x... i read for the fast FFT it doesnt matter.. better do it on the CPU. or where ever the data is ... visual studio update has some texture compression stuff .. and i read if its 64k i think.. its fits in one warp.... i might try Cudaify.. .. i was doing almost realtime mandelbrot at full resolution...with an opencl demo.. so im planning ot use some 1 bit.. 2 bit .. paletted bitmaps..if i have to pass them around.

i would love if openCL finds the processor in skylake, and if the is an AMD or nvidia ,,the we could use them both on newer machines.. skylake GPU on the chip is almost good enough,, not quite ... so alot of systems will have two powerful gpus. i cant seem to query my intel one..but i think its too old..

im still thinkingo of the FFT solv.. there are alot of tricks.. here are a few more papers.. if its just for pressure solv in air...the presence of one small 2d character, would not ruin it..just not be perfectly accurate? the FFT region might just be the sky ... then using some other grids in chimera style....

i guess its pretty easy to try them... for a side scoller game.. the boundaries are often good enough..

http://ms.mcmaster.ca/~bprotas/MATH745c/spectr_01.pdf

http://ufrmeca.univ-lyon1.fr/~buffat/FDMIF/articleIBM.pdf

.https://people.eecs.berkeley.edu/~demmel/cs267_Spr05/Lectur es/Lecture26/Lecture_26_AMR_part2_PhilColella_05.pdf

On Mon, Sep 12, 2016 at 5:06 AM, Gonçalo Amador notifications@github.com wrote:

Used GPU before with CUDA and OpenCL the times in the thesis actually can be way faster.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/noSoulApophis/incremental-fluids-java/issues/2#issuecomment-246204629, or mute the thread https://github.com/notifications/unsubscribe-auth/AP_LmaSlxuQrAmU6MHvW0vp0Sgat3H8Fks5qpG1RgaJpZM4J56kt .

[damian-666]

just ignore these if they make no sense i am still learning. and i am only interested in 2d physics for my project.​

BTW .. if you havent seen this .. its MIT GPU + java + processing... https://github.com/diwi/PixelFlow im going to try it but it might suit you more than me.. http://thomasdiewald.com/blog/

i sent this because you didnt like the FFT method with IB .. and i have a hard time also to believe, since its is so much faster on large grids than others which we see normally used .. the MG.. the preconditioners.. the CG.. the interative solvers... but i find here a nice paper again..​ and it makes sense. the idea is the BC must be periodic at left and right and top and bottom not sure.. but we can hide that in the margin of the viewport. if the flow leaves the viewport.. we can force it in the Y direction so that it doesnt resemble what left... for my top and bottom i would put a margin.. and have a "jet stream" or fast wind with low shear. for my bottom i would have an ocean flow or some simply decribed sine shaped hills ( easiest for the solver)

One can use the liquid surfaces as an Immersed boundaries as well. use a level set represent it..advence that level set and use a rough grid to the bottom or maybe the shallow equations.. or... solve a separate grid for the water..

A ​nd here is the central idea.. in the 2015​

​IMERSPECMETHODOLOGYFORTWO-PHASEFLOW

​The Fourier pseudospectral method, is a methodology of high rate of numerical convergence and has high accuracy and low computational cost when compared with other methodologies of high accuracy due to the use of the algorithm namedFastFourierTransform(FFT).ItisalsoobservedthatfortheNavier-Stokesequations,consideringincompressible 1The term accuracy is used in this report as the difference between the numerical solution and real solution.

flow in spectral space, the projection method can be used to eliminate the pressure gradient term, leading to uncoupling of the pressure-velocity fields without the need to solve the Poisson equation for pressure. The major limitation of this methodology is in the boundary conditions, which must be periodical in order to apply the Fourier spectral method, (Canuto et al., 2007). However, Mariano et al. (2010) overcame this limitation by coupling the Fourier pseudospectral method with the immersed boundary method, which allowed to simulate non-periodic problems, using pseudospectral method. The search ​ which allowed to simulate non-periodic problems, using pseudospectral method. The search for more accurate and low computational cost for the simulation of two-phase flow methods is of great interest for industry, since they require information from the flows with higher level of accuracy. This justifies financial support and technical cooperation of Petrobras in the present work. Thus, the purpose of the present paper is to show a proposition of a new mathematical and computational methodology to solve two-phase flows gas-liquid, liquid-liquid or gas-gas problems which provides at the same time, high computational efficiency and high accuracy. For this we used Fourier pseudospectral method (FSM) coupled with immersed boundary (IB) and with the Front-Tracking method (FTM). Results obtained by analysis through numerical experiment of rising bubbles show that the proposed method can be considered validated and very promising.

2. MATHEMATICALMODELINGFORTWOPHASEFLOWS The present work is based on the merging FSM, IB and FTM method. As shown in Fig. 1, the fluid as a whole is representedby Ω = Ω1∪Ω2 domain. Theinterface Γ iscalledLagrangian,whichcanmoveanddeform. Thesubdomains Ω1 and Ω2 represents the outside and inside of the interface Γ, respectively. Variables with uppercase letter are related to Lagrangian domain and variables with lowercase letter are related to Eulerian field.