Water-floater(solid structure) interaction

[MOHAMED EL AMINE LAKHNICHI]

I need to simulate the interaction between water(incompressible) with a solid structre(simple floater) in order to evaluate the energy transport between the two systems.I've just known about Pysph so i don't know from where I should start. If you can help please !

[Marina Kauffeldt]

There examples for your problem in the documentation https://pysph.readthedocs.io/en/latest/examples/sphere_in_vessel.html

and there are more in the git repository https://github.com/pypr/pysph/tree/master/pysph/examples/rigid_body

good luck :)

[A Dinesh]

Hi Mohamed,

I have implemented rigid fluid coupling in PySPH. Your problem specifies you want to simulate interaction between a structure and fluid, assuming structure here to be rigid I would like to continue. I would like to answer this question in steps and use it as a future documentation for rigid fluid coupling in PySPH. So here we go.

I assume you know some python (no other prerequisites for now to do the coupling problem). Since you don't know where to start, I would like you to first install the software. Here are the installation instructions. After installation please try to run the examples (This is check

if the software is installed correctly. I would ask to not run tests as some tests might fail, and it doesn't stop us from what we are going to doing).

Before jumping right away to the coupling problem, lets first solve some single phase problems (I mean only fluid). Some concrete examples to start with would be at hydrostatic tank, dam break 2d and dam break 3d. Running these examples will give you some exposure to the PySPH sofware and also some hand on python. There are two levels of doing this (actually three)

1. Running directly using `pysph` command
2. Running using the python file
3. Running by writing our own example
   

For now I will explain the first two ways.

## Running directly using the `pysph` command

Assuming you are on Linux based OS (I don't know how these commands work on  Windows), if you open a `terminal` and having `pysph` installed, you can do,

pysph run

this will display all the example `pysph` can run, it will be something like following,

|  LAPTOP LAPTOP=> pysph run
1. flow_past_cylinder_2d
   A wind-tunnel flow past a cylinders.  (20 minutes)
2. two_blocks
   Two square blocks of water colliding with each other. (20 seconds)
3. hydrostatic_tank
   Hydrostatic tank example. (2 minutes)
4. couette
   Couette flow using the transport velocity formulation (30 seconds).
5. rayleigh_taylor
   Rayleigh-Taylor instability problem. (16 hours)
6. cavity
   Lid driven cavity using the Transport Velocity formulation. (10 minutes)
   .
   .
   .
   .
   .
   .
34. solid_mech.impact3d
   High-velocity impact of an Steel projectile on an Aluminium plate"""
35. solid_mech.impact
   High-velocity impact of an Steel projectile on an Aluminium plate"""
36. solid_mech.post_process_oscillating_plate
   ort numpy as np
37. solid_mech.oscillating_plate
   ort numpy as np
38. solid_mech.rings
   Colliding Elastic Balls. (10 minutes)
39. solid_mech.oscillating_plate_gtvf
   ort numpy as np
40. solid_mech.taylor_bar
   Taylor bar example with SPH. (5 minutes)

You may not see the same number of examples. I am using a different branch (you should know about `git` to understand this, if you don't, you don't have to worry) on my machine. On my machine the `dam break 2d` example is number 13,

by typing `13` and enter, the example will run. It would look like

Enter example number you wish to run: 13
Enter additional arguments (leave blank to skip):
--------------------------------------------------------------------------------
Running example pysph.examples.dam_break_2d.

Information for example: pysph.examples.dam_break_2d
Two-dimensional dam break over a dry bed.  (30 minutes)

The case is described in "State of the art classical SPH for free surface
flows", Moncho Gomez-Gesteira, Benedict D Rogers, Robert A, Dalrymple and Alex
J.C Crespo, Journal of Hydraulic Research, Vol 48, Extra Issue (2010), pp
6-27. DOI:10.1080/00221686.2010.9641242
Using h = 0.039000
2D dam break with 2211 fluid, 1624 boundary particles
Generating output in /home/dinesh/dam_break_2d_output
----------------------------------------------------------------------
No of particles:
  fluid: 2211
  boundary: 1624
  Total: 3835
----------------------------------------------------------------------
Setup took: 13.14449 secs
 25%|█████               | 3.4kit | 6.1e-01s [00:23.1<01:10.8 | 0.007s/it]

Assignment 1: Please run this example on your system and print the following.

1. Output folder name
2. Total time taken to run the simulation
3. Post the image of the output of your simulation. (This can be done by taking screen shot of the `pysph viewer`)

SImilarly you can run the other examples.

Summary:

I have asked to do the following,

1. Install the software
2. Run the examples with `pysph run` command. (This is the simplest of above three ways of running. This shouldn't be hard)
3. Do the Assignment 1.
   
4. This should not take more that an hour. 30 minutes for installation + 30 for running the simulation and doing the assignment 1. 

Please do this and reply. If you found any difficulties in doing any of the above please ask here as soon as possible. This would be good exercise to see how I can taper the tutorial for wider audience.

Regards,

Dinesh.

[MOHAMED EL AMINE LAKHNICHI]

Hi Dinesh,

Thanks a lot for the initiative. Anyway, I'm working on windows (10) and tu run Pysph I use Windows prompt; I use that path c:\programdata\anaconda3\Scripts where I installed

 python-anaconda, pysph as well. Then I run the example taylor_green. Here's the result I found

After running the examples, nothing was showed up, I tried then to use pysph view and it doesn't work

Mohamed,

Thanks

PS sorry for the late reply

[MOHAMED EL AMINE LAKHNICHI]

[MOHAMED EL AMINE LAKHNICHI]

Hi,

I used the indicated instructions in the website about sph, and it works:

[MOHAMED EL AMINE LAKHNICHI]

Le dimanche 19 janvier 2020 16:11:40 UTC+1, MOHAMED EL AMINE LAKHNICHI a écrit :

[MOHAMED EL AMINE LAKHNICHI]

Hi Dinesh,

You asked me for the dam tank 2d and 3d, but the 3d requires 14 hours

--
You received this message because you are subscribed to the Google Groups "pysph-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/50569a38-c764-442c-8c17-e1078c2a126f%40googlegroups.com.

[A Dinesh]

For now please run the other two examples. You don't have to run them for long time (if the simulation is too long), after running for some permissible time, stop it and post the images.  Please post the name of the output directory.

On Sunday, February 2, 2020 at 5:51:25 PM UTC+5:30, MOHAMED EL AMINE LAKHNICHI wrote:

Hi Dinesh,

You asked me for the dam tank 2d and 3d, but the 3d requires 14 hours

Le dim. 2 févr. 2020 à 13:11, MOHAMED EL AMINE LAKHNICHI <lakhnichi.mohamedelamine.mp...@gmail.com> a écrit :

Le dimanche 19 janvier 2020 16:11:40 UTC+1, MOHAMED EL AMINE LAKHNICHI a écrit :

I need to simulate the interaction between water(incompressible) with a solid structre(simple floater) in order to evaluate the energy transport between the two systems.I've just known about Pysph so i don't know from where I should start. If you can help please !

--
You received this message because you are subscribed to the Google Groups "pysph-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph...@googlegroups.com.

[MOHAMED EL AMINE LAKHNICHI]

How could I know that output directory, sorry if I didn't understand it well!

Le dimanche 19 janvier 2020 16:11:40 UTC+1, MOHAMED EL AMINE LAKHNICHI a écrit :

[MOHAMED EL AMINE LAKHNICHI]

I run the dam tank -2d version- here are the results:

Still don't know about the output directory, but i would like to know if it is possible to play the simulation automatically ?

Thanks,

Le dimanche 19 janvier 2020 16:11:40 UTC+1, MOHAMED EL AMINE LAKHNICHI a écrit :

[A Dinesh]

Output directory is something, where your data is dumped at an intervals of time. Here the output folder is `dam_break_2d_output`.

[A Dinesh]

Next, we will run some examples using python command. If you have installed `PySPH` using pip, please download the `PySPH` repository on your local machine.

The package is here. Or if you have installed it through repository using the following command,

python setup.py install

then change into the directory. The following illustration is based on Linux OS. First change to the directory of `PySPH`, then goto examples directory. This is done as

cd pysph

If you do `ls` here you should see the following,

$LAPTOP LAPTOP=> ls
appveyor.yml  docs         old_examples    requirements-test.txt  shippable.yml
build         LICENSE.txt  pysph           requirements.txt       starcluster
CHANGES.rst   Makefile     PySPH.egg-info  setup.cfg              tox.ini
docker        MANIFEST.in  README.rst      setup.py

Now goto directory `pysph/examples`.

$ LAPTOP LAPTOP=> cd pysph/examples

$ LAPTOP LAPTOP=> ls
cavity.py                     gas_dynamics           run.py
couette.py                    ghia_cavity_data.py    run.pyc
cube.py                       hydrostatic_tank.py    solid_mech
dam_break_2d.py               __init__.py            spheric
dam_break_3d.py               __init__.pyc           sphysics
db_exp_data.py                lattice_cylinders.py   surface_tension
_db_geometry.py               periodic_cylinders.py  taylor_green_output
dem                           poiseuille.py          taylor_green.py
elliptical_drop_no_scheme.py  __pycache__            tests
elliptical_drop.py            rayleigh_taylor.py     trivial_inlet_outlet.py
elliptical_drop_simple.py     rigid_body             two_blocks.py
flow_past_cylinder_2d.py      rigid_fluid

Now try running some examples with `python` command

python dam_break_2d.py

There is no need to understand what is inside the code for now, we will go in detail, how to write our own simulations and physics.

Try running the same example with few command line options `PySPH` provides us.

python dam_break_2d.py --openmp

or do

python dam_break_2d.py --opencl

if you have GPU card on your machine.

PySPH has quite a few implementations to simulate fluid flow, such as `WCSPH`, `TVF`, `ISPH`, `IISPH` etc. Try running the same example with three different schemes

python dam_break_2d.py --openmp --scheme wcsph

Check all the command line arguments `PySPH` provides to run a simuation.

To summarize,

1. Download (or clone (for people who know git))  the repository.
2. Go to the examples directory.
3. Run a few examples using `python` command.
   
4. Try few different schemes by changing the command line argument
   
5. Try three different examples with above arguments.
6. Try running with `openmp` and `opencl` command line arguments

Assigment:

1. Run three examples with python command
2. Run with `openmp` and `opencl` (`opencl` only if you have a GPU card)
3. Play around with the command line arguments (hint: you can see all the arguments with `python example_name.py --help`

Good luck.

[MOHAMED EL AMINE LAKHNICHI]

Hi,

I faced a problem after trying to run the example using python, I switched to another; the execution was accomplished but once I tried use the '--openmp' command it raises an error

Le dimanche 19 janvier 2020 16:11:40 UTC+1, MOHAMED EL AMINE LAKHNICHI a écrit :

[Enrique del Castillo]

Hi Dinesh,

Thank you for taking the time to organize this tutorial. I have also been following the steps and was able to complete the first task with the 2d dam.

Next in the second task I did:

$ git clone https://github.com/pypr/pysph.git

$ cd pysph 

$ cd pysph/examples

$ python dam_break_2d.py

and seem to be running into the same error that Mohamed is:

I next tried 

$ python elliptical_drop.py

getting the following output:

When I navigate to the elliptical_drop_output folder I see this:

I am able to open comparison.png, but I cannot open any of the other files, using for example,

$ pysph view elliptical_drop_300.hdf5

or even using like in task 1 just,

$ psych view elliptical_drop_output

How should we be visualizing this output?

I am also a little bit confused by the fact that it seems we have created two different pysph folders, one created in part 2 with the address:

/Users/EnriquedelCastillo/pysph/

and one from part 1:

/tmp/pysph/

Yet here the elliptical_drop_output folder was in the address from part 1

/tmp/pysph/pysph/examples/elliptical_drop_output

whereas 

/Users/EnriquedelCastillo/pysph/pysph/examples/elliptical_drop_output

does not exist. I was wondering why the elliptical_drop_output file is going to this address, or if maybe I have done something wrong in the installation steps perhaps?

Thank you!

[Enrique del Castillo]

For some reason the pictures did not load in my original post. Here they are in order.

[A Dinesh]

I guess you are might have installed PySPH using pip. Since pip doesn't provide bleeding edge version of PySPH, I would ask you to install through repository.

1. Uninstall `PySPH` if you have installed it through `pip`.
2. Install from the repository itself.
   

To install from the repository, first install the dependencies of `PySPH` this can be done with

pip instal -r requirements.txt

After that install the software by

python setup.py develop

or

python setup.py install

prefer `develop` over `install`, as we will see why in the following posts.

After doing the above please try to complete the assignment.

[MOHAMED EL AMINE LAKHNICHI]

The directory is the one from https://github.com/pypr/pysph.git ?

[Enrique del Castillo]

Hi Dinesh,

Yes I had installed it with pip and then followed the installation steps using git from the github repository itself, but I had not uninstalled the pip PySPH installation before doing that. 

I have now uninstalled the pip PySPH and deleted the pysph folder from the git cloning, and redone the installation from the repository using git and 

pip install -r requirements.txt

as well as

python setup.py develop

as you suggested, although I had already installed the “Core dependencies” in the PySPH documentation as well as those in the “Installing the dependencies on Mac OS X” and “Using Anaconda” sections of the document.

Now however, when I try running $python dam_break_2d.py

I receive a segmentation fault error after 38% progress is reached:

(Image 1)

I also tried $python setup.py install, but achieved the same result.

For $python dam_break_2d.py —openmp

I get the subsequent error:

(Image 2)

Thanks for the help.

[A Dinesh]

The idea of running the examples with `--openmp` flag is to use all the cores one system provides. Basically with `--openmp` it will be a parallel code. (Fix instructions will be posted soon).

The example on master branch is failing with the default scheme (which is PySPH). Try running the example with different scheme (such as tvf, or aha or edac).

And also more two examples so that you will get some practice, with different command line arguments.

[Prabhu Ramachandran]

Hi,

Here are a few things that may help:

1. I think the dam break 2d problem in master is broken for the WCSPH case as for that scheme the particle packing needs to be hexagonal ideally.  So for now don't worry about this.  There are things in git master that need fixing that we have not done a good job of fixing in time.  Sorry about that.

2. On MacOS (which I also use), setting up openmp can be a bit tricky.  How have you installed gcc?  Is it via conda or via brew?  I prefer using brew.  Make sure there aren't two different and incompatible gcc installations, for example conda can also install gcc and I have not tested with that configuration.  It is possible that conda's gcc and openmp configuration is in the way.  The pysph instructions should hopefully cover the openmp installation on MacOS via brew.

Regards,

Prabhu

p.s. You don't need to send screenshots of your terminal and it may be more reliable to just send it as text in the email.

--

You received this message because you are subscribed to the Google Groups "pysph-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/375861e6-d280-4d2c-b860-92a3a8f55561%40googlegroups.com.

-- 
Prabhu Ramachandran                 http://www.aero.iitb.ac.in/~prabhu

[Enrique del Castillo]

Hi,

I am unfortunately still unable to use openmp, and get the same error as before. I did already have gcc installed, and I had followed the steps in the pysph instructions using brew. $ brew install gcc $ export CC=gcc-9  $ export CXX=g++-9. Is this all that was needed for the gcc installation for MacOS? If I check my gcc version by  typing g++ and pressing tab twice, I see g++-9 displayed.

With regards to the dam_break_2d, I tried out some of the different schemes and they work with no openmp for example:

$ python dam_break_2d.py --no-openmp --scheme aha

$ python dam_break_2d.py --no-openmp --scheme edac

$ python dam_break_2d.py --no-openmp --scheme sisph

$ python dam_break_2d.py --no-openmp --scheme gtvf

I have been able to run and visualize some other examples such as the poiseuille.py, couette.py, and elliptical_drop.py without any trouble. I experimented using different arguments such as:

$ python poiseuille.py --no-openmp --alpha .5

$ python poiseuille.py --no-openmp --tdamp 2

$ python couette.py --no-openmp --nnps sh

$ python couette.py --no-openmp --zoltan-weights 

$ python elliptical_drop.py --iisph --tensile-correction

among others, encountering no issues. 

Is it possible to proceed onto the next steps without openmp (hopefully temporarily, while I continue to try to fix it)?

Thanks!

To unsubscribe from this group and stop receiving emails from it, send an email to pysph...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/375861e6-d280-4d2c-b860-92a3a8f55561%40googlegroups.com.

[A Dinesh]

I think there no option like `no-openmp`. I will post the next task by tonight.

[Prabhu Ramachandran]

On 2/5/20 5:14 AM, Enrique del Castillo wrote:

Hi,

I am unfortunately still unable to use openmp, and get the same error as before. I did already have gcc installed, and I had followed the steps in the pysph instructions using brew. $ brew install gcc $ export CC=gcc-9  $ export CXX=g++-9. Is this all that was needed for the gcc installation for MacOS? If I check my gcc version by  typing g++ and pressing tab twice, I see g++-9 displayed.

With regards to the dam_break_2d, I tried out some of the different schemes and they work with no openmp for example:

$ python dam_break_2d.py --no-openmp --scheme aha

$ python dam_break_2d.py --no-openmp --scheme edac

$ python dam_break_2d.py --no-openmp --scheme sisph

$ python dam_break_2d.py --no-openmp --scheme gtvf

Do you have any openmp libraries (or gcc) inside your conda environment?  If so this could be a problem as it is possible that this library is being linked.  One option is to not export CXX and try (that should pick up any local GCC).  For example try this:

$ conda list | grep openmp

and if this returns something or the other then it is possible that something is not linked correctly.  If this is the case, perhaps forget the brew installation and try with the conda version of gcc and see if that works. 

If you really want to get to the bottom of this, you can look for the extension module that is compiled (see your error message), for example on my system I have this:

$ otool -L ~/.pysph/source/py3.6-macosx-10.6-x86_64/m_189951b83f420e1adcb301ba3c31024c.cpython-36m-darwin.so

[...]

   /usr/local/opt/gcc/lib/gcc/8/libgomp.1.dylib (compatibility version 2.0.0, current version 2.0.0)
    /usr/local/lib/gcc/8/libgcc_s.1.dylib (compatibility version 1.0.0, current version 1.0.0)

As you can see, it links to the libgomp from brew and not anywhere else.

Finally, if you do not use the --openmp argument, it will NOT use openmp and everything ought to work without using OpenMP.

cheers,

Prabhu

[A Dinesh]

Till now we ran our simulations using `PySPH run` command and by using `schemes` though `python example.py --scheme gtvf` command. Now we will look at the third way of running the simulation i.e., without using schemes.

In the examples folder look at `elliptical_drop_no_scheme.py` example, run it.

Do the tutorial ipython notebook (Explains about writing simulations without using scheme). This gives a nice overview about what methods we need overload in order to run a simulation. In the upcoming part we will look at each individual method.

I will explain creating particle array (fluid, or solid, etc), what are groups in equations and what is a solver etc.

[A Dinesh]

Looks like these (link) tutorials cover basics of PySPH. Before going forward, I want you guys to complete those tutorials. I will explain what not covered in those ipynotebooks. Finally we will implement discrete element method on our own in `PySPH` and run hopper flow. 

I will explain how equations work (this is not covered in notebooks)  in the next part. I have a simple assignment to start with. Consider a n-body problem ( I hope you know how a n body problem works, if not please ask, I will try to post more resources). Rosetta code has python implementation here. Look at the `_resolveCollisions` function and try to write your own `resolveCollision` (It is not much different) and post here. This will be a good start to understand the upcoming next part. Next part will come tonight after some one completes the above tasks.

On Wednesday, February 5, 2020 at 5:14:51 AM UTC+5:30, Enrique del Castillo wrote:

[A Dinesh]

A little contest for the above task. Create a particle array (You will be able to understand this after completion of ipython notebooks) with 10 particles, and name it as fluid. Assuming each particle influences every other particle (N-body problem), then what would be the force acting on each particle if force follows the following equation,

f_ij = k_i * m_i * m_j / (rij_2)

f_ij is force on i due to j. Where k_i is constant of each particle, m_i and m_j are masses and (r_ij) is the distance between i and j particles. Try to write this in a pure python for loop.

Assingment:

Write a simple interaction loop in pure python, and find the force on all particles.

[Enrique del Castillo]

Hi,

Sorry for the delay in replying. I've been having some issues again running some of the files in the tutorial. For example after I checked my solutions with the ed.py and the ed_no_scheme.py, I saved these files in the same examples folder and tried running them using

$ python3 ed.py

for example, and received the following error and warnings which are confusing me a bit.

Generating output in /Users/EnriquedelCastillo/pysph/pysph/examples/ed_output

********************************************************************************

ERROR

running build_ext

cythoning /Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.pyx to /Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.cpp

building 'm_eca5403f6fd2d9e0c9b4331cfa431069' extension

/usr/bin/clang -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/Users/EnriquedelCastillo/.conda/envs/py3/include -arch x86_64 -I/Users/EnriquedelCastillo/.conda/envs/py3/include -arch x86_64 -I/Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/numpy/core/include -I/Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/PySPH-1.0b1.dev0-py3.5-macosx-10.6-x86_64.egg/pysph/base -I/usr/local/include -I/Users/EnriquedelCastillo/.conda/envs/py3/include/python3.5m -c /Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.cpp -o /Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/build/temp.macosx-10.6-x86_64-3.5/Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.o

/Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/Cython/Compiler/Main.py:369: FutureWarning: Cython directive 'language_level' not set, using 2 for now (Py2). This will change in a later release! File: /Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.pyx

  tree = Parsing.p_module(s, pxd, full_module_name)

warning: /Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/PySPH-1.0b1.dev0-py3.5-macosx-10.6-x86_64.egg/pysph/base/nnps_base.pxd:25:13: 'INT_MAX' redeclared 

warning: include path for stdlibc++ headers not found; pass '-stdlib=libc++' on the command line to use the libc++ standard library instead [-Wstdlibcxx-not-found]

In file included from /Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.cpp:613:

In file included from /Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/numpy/core/include/numpy/arrayobject.h:4:

In file included from /Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/numpy/core/include/numpy/ndarrayobject.h:12:

In file included from /Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/numpy/core/include/numpy/ndarraytypes.h:1822:

/Users/EnriquedelCastillo/.conda/envs/py3/lib/python3.5/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: "Using deprecated NumPy API, disable it with "          "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-W#warnings]

#warning "Using deprecated NumPy API, disable it with " \

 ^

/Users/EnriquedelCastillo/.pysph/source/py3.5-macosx-10.6-x86_64/m_eca5403f6fd2d9e0c9b4331cfa431069.cpp:615:10: fatal error: 'ios' file not found

#include "ios"

         ^~~~~

2 warnings and 1 error generated.

Is it a problem with my cython or gcc again? Wierdly, when I try running any of the example scripts such as couette.py etc., no such error appears and they run fine.

Also, might there be a solution for the advanced exercise in tutorial2?

Thank you!

[Prabhu Ramachandran]

On Thursday 06 February 2020 07:12 AM, A Dinesh wrote:

Looks like these (link) tutorials cover basics of PySPH. Before going forward, I want you guys to complete those tutorials. I will explain what not covered in those ipynotebooks. Finally we will implement discrete element method on our own in `PySPH` and run hopper flow. 

I will explain how equations work (this is not covered in notebooks)  in the next part. I have a simple assignment to start with. Consider a n-body problem ( I hope you know how a n body problem works, if not please ask, I will try to post more resources). Rosetta code has python implementation here. Look at the `_resolveCollisions` function and try to write your own `resolveCollision` (It is not much different) and post here. This will be a good start to understand the upcoming next part. Next part will come tonight after some one completes the above tasks.

Please note that true n-body problems cannot be solved with pysph currently.  What you can do is use a brute-force N^2 algorithm with nearest neighbors suitably setup. 

Prabhu

[Prabhu Ramachandran]

On Sunday 09 February 2020 12:34 PM, Enrique del Castillo wrote:

Hi,

Sorry for the delay in replying. I've been having some issues again running some of the files in the tutorial. For example after I checked my solutions with the ed.py and the ed_no_scheme.py, I saved these files in the same examples folder and tried running them using

$ python3 ed.py

for example, and received the following error and warnings which are confusing me a bit.

You seem to be using clang and for some reason clang is not able to find ios?  I suspect you forgot to set the CXX env variable.  The other examples may work because they were compiled earlier and it is now using the cached extension module.

Prabhu

[A Dinesh]

You try them out. If you have any questions, ask here. I will try to answer.

[A Dinesh]

Hi Enrique,

Does this folder provide solutions? Please check.

[Enrique del Castillo]

Hi,

Thank you Prabhu, that was exactly the issue! I have now fixed it. 

Dinesh - 

Here is my script for the interaction loop:

from __future__ import print_function

from pysph.base.particle_array import ParticleArray

from matplotlib import pyplot as plt

pa = ParticleArray(name='fluid', x=[1,2,3,4,5,6,7,8,9,0], y=[.3,5,1,8,5,3,1,3,3,5],m=[0.1,2,1,3,4,5,6,7,2,4],k=[0.1,.3,.6,.35,.43,.5,.6,.7,.21,.9] )

plt.scatter(pa.x, pa.y, marker='.') 

plt.xlabel('x')

plt.ylabel('y')

force_list = [] #list of the total forces for each particle i (Tf_i) to be printed

for i in range(len(pa.m)):

    Tf_i = 0 #total force on each particle i

    for j in range(len(pa.m)):

        #exclude case where i=j, cannot have force on particle produced by itself

        if j!=i:

            k_i = pa.k[i]

            m_i = pa.m[i]

            m_j = pa.m[j]

            rij = ((pa.x[j]-pa.x[i])**2)+((pa.y[j]-pa.y[i])**2)**.5

            rij_2 = rij**2 # I assumed underscore 2 means squared.

            f_ij = k_i*m_i*m_j/(rij_2)

            Tf_i = Tf_i + f_ij

    force_list.append(Tf_i)

print(force_list)

I had seen the solution folder within the tutorial which you put in the link, but I meant to ask about the dam_break.py. "advanced example" at the bottom of tutorial 2, I don't think it is in the solution folder. In any case it is Ok, I am ready to continue with the Hopper flow task.

Thank you.

[A Dinesh]

On Monday, February 10, 2020 at 3:02:06 AM UTC+5:30, Enrique del Castillo wrote:

Hi,

Thank you Prabhu, that was exactly the issue! I have now fixed it. 

Dinesh - 

Here is my script for the interaction loop:

from __future__ import print_function

from pysph.base.particle_array import ParticleArray

from matplotlib import pyplot as plt

pa = ParticleArray(name='fluid', x=[1,2,3,4,5,6,7,8,9,0], y=[.3,5,1,8,5,3,1,3,3,5],m=[0.1,2,1,3,4,5,6,7,2,4],k=[0.1,.3,.6,.35,.43,.5,.6,.7,.21,.9] )

plt.scatter(pa.x, pa.y, marker='.') 

plt.xlabel('x')

plt.ylabel('y')

force_list = [] #list of the total forces for each particle i (Tf_i) to be printed

Nice. I just have a few suggestions to the above code. Rather than appending the force to a new list i.e., force_list.append(Tf_i), why don't  you create a property called `Tf_i` (you should be able to do this if you have done particle array tutorial) to the particle array itself, and before starting the loop make the forces equate to zero,

then add the force to the respective particle. Why I am asking this is, we follow similar philosophy in PySPH.

 

for i in range(len(pa.m)):

    Tf_i = 0 #total force on each particle i

    for j in range(len(pa.m)):

        #exclude case where i=j, cannot have force on particle produced by itself

        if j!=i:

            k_i = pa.k[i]

            m_i = pa.m[i]

            m_j = pa.m[j]

            rij = ((pa.x[j]-pa.x[i])**2)+((pa.y[j]-pa.y[i])**2)**.5

            rij_2 = rij**2 # I assumed underscore 2 means squared.

            f_ij = k_i*m_i*m_j/(rij_2)

            Tf_i = Tf_i + f_ij

    force_list.append(Tf_i)

print(force_list)

I had seen the solution folder within the tutorial which you put in the link, but I meant to ask about the dam_break.py. "advanced example" at the bottom of tutorial 2, I don't think it is in the solution folder. In any case it is Ok, I am ready to continue with the

Ahh, okay. If the solutions are not provided, please try it yourself and post any questions here. I will try to answer them.

Hopper flow task.

I am touching on all the required topics needed to know before we go on implementing the algorithm our self. This may take more one week. Depends on how fast we go.

[Enrique del Castillo]

Hi,

 I have fixed the script as suggested. I noted that by leaving the array as Tf_i = [] , this produces values of 0, and makes the array the length of the other particlearray property arrays. 

from __future__ import print_function

from pysph.base.particle_array import ParticleArray

from matplotlib import pyplot as plt

pa = ParticleArray(name='fluid', x=[1,2,3,4,5,6,7,8,9,0], y=[.3,5,1,8,5,3,1,3,3,5],m=[0.1,2,1,3,4,5,6,7,2,4],k=[0.1,.3,.6,.35,.43,.5,.6,.7,.21,.9], Tf_i = [])

plt.scatter(pa.x, pa.y, marker='.') 

plt.xlabel('x')

plt.ylabel('y')

for i in range(len(pa.m)):

    Tf_i_temp = 0 #total force on each particle i

    for j in range(len(pa.m)):

        #exclude case where i=j, cannot have force on particle produced by itself

        if j!=i:

            k_i = pa.k[i]

            m_i = pa.m[i]

            m_j = pa.m[j]

            rij = ((pa.x[j]-pa.x[i])**2)+((pa.y[j]-pa.y[i])**2)**.5

            rij_2 = rij**2 

            f_ij = k_i*m_i*m_j/(rij_2)

            Tf_i_temp = Tf_i_temp + f_ij

    pa.Tf_i[i] = Tf_i_temp

I was able to finish that dam_break.py example. I based my script below on ed.py and I took the boundary geometry  from the actually dam_break_2d.py example.

import numpy as np

from pysph.base.utils import get_particle_array

from pysph.solver.application import Application

from pysph.sph.scheme import WCSPHScheme

from pysph.tools.geometry import get_2d_tank

rho = 1.0

g = 9.81

container_height = 4.0

container_width = 4.0

nboundary_layers = 2

dx = 0.025

class DamBreak(Application):

    def create_particles(self):

        x, y = np.mgrid[0:1:dx, 0:2:dx]

        xt, yt = get_2d_tank(dx=dx, length=container_width,

                             height=container_height, base_center=[2, 0],

                             num_layers=nboundary_layers)

        h = 1.3*dx

        m = rho*dx*dx

        fluid = get_particle_array(

            name='fluid', x=x, y=y, rho=rho,

            m=m, h=h

        )

        boundary = get_particle_array(

            name='boundary', x=xt, y=yt, rho=rho,

            m=m, h=h

        )

        self.scheme.setup_properties([fluid, boundary])

        return [fluid, boundary]

    def create_scheme(self):

        s = WCSPHScheme(

            ['fluid'], ['boundary'], dim=2, rho0=1.0, c0=1400,

            h0=1.3*0.025, hdx=1.3, gamma=7.0, gy=-g, alpha=0.1, beta=0.0

        )

        dt = 5e-6

        tf = 2.5

        s.configure_solver(

            dt=dt, tf=tf,

        )

        return s

if __name__ == '__main__':

    app = DamBreak(fname='dam_break')

    app.run()

The behavior of the fluid is quite different from the dam_break_2d.py example however, and looks physically impossible especially along the left boundary. I put some pictures attached. Did I make some mistake, or is it related to the scheme used? (Vmagnitude is plotted).

Thanks!

[A Dinesh]

On Monday, February 10, 2020 at 1:27:28 PM UTC+5:30, Enrique del Castillo wrote:

Hi,

 I have fixed the script as suggested. I noted that by leaving the array as Tf_i = [] , this produces values of 0, and makes the array the length of the other particlearray property arrays. 

from __future__ import print_function

from pysph.base.particle_array import ParticleArray

from matplotlib import pyplot as plt

pa = ParticleArray(name='fluid', x=[1,2,3,4,5,6,7,8,9,0], y=[.3,5,1,8,5,3,1,3,3,5],m=[0.1,2,1,3,4,5,6,7,2,4],k=[0.1,.3,.6,.35,.43,.5,.6,.7,.21,.9], Tf_i = [])

Or you can simply do `tf_i=0`, this would initialize all the values with 0 (with array of length of other particle array properties). There is one more option

`pa.add_property('tf_i`)`

If you want the property to be of type `int`, then you do

`pa.add_property('tf_i`, type='int')`

this will be useful when implementing some different physics ("discrete element method").

I have an interesting question here. Say I have a particle array as

x = np.linspace(0., 1., 10)
pa = get_particle_array(x=x, h=1., m=1.) # this function is in pysph.base.utils


# what do you think the following will have?
pa.m  # will this be an array or a float
pa.h  # same question here

# Question 2: similarly
pa = get_particle_array(x=x, h=1., m=[1.]) # this function is in pysph.base.utils

# what do you think the following will have?
pa.m  # will this be an array or a float Question (2a)
pa.h  # same question here Question (2b)

# Question 3: similarly
pa = get_particle_array(x=x, h=1., m=[1., 2.]) # this function is in pysph.base.utils


# what do you think the following will have?
pa.m  # will this be an array or a float Question (3a)
pa.h  # same question here Question (3b)


Sometime you might add the particle property and intialize it as follows,

x = np.linspace(0., 1., 10)
pa = get_particle_array(x=x, h=1., m=1.)
pa.add_property('tf_i') # We now know that this is initialized with 10 length array (cyarray) with 0 values
# but I want to change the value to 10. now.


# Question 4
pa.tf_i = 10. # What happens here (question 4a)

pa.tf_i[:] = 10. # What happens here (question 4b)

plt.scatter(pa.x, pa.y, marker='.') 

plt.xlabel('x')

plt.ylabel('y')

Assuming that the particle property total force ('tf_i`) is not zero, then we would do,

for i in range(len(pa.m)):

    pa.tf_i[i] = 0

Here we don't need to run this single for loop. But supposing we are in a simulation, to compute the force at the current time step, we should always make it zero and start over (I hope you know this).

Then we will compute the force as you wrote.

for i in range(len(pa.m)):

I have a few more suggestions. Idiomatic PySPH code would be like

    # we don't need this line

    # Tf_i_temp = 0 # total force on each particle i

    for j in range(len(pa.m)):

        #exclude case where i=j, cannot have force on particle produced by itself

        if j!=i:

            k_i = pa.k[i]

            m_i = pa.m[i]

            m_j = pa.m[j]

            rij = ((pa.x[j]-pa.x[i])**2)+((pa.y[j]-pa.y[i])**2)**.5

            rij_2 = rij**2 

            f_ij = k_i*m_i*m_j/(rij_2)

  DIrectly add the force to the particle property.

            pa.Tf_i[i] _+= f_ij

    # pa.Tf_i[i] = Tf_i_temp

And also remove the above line.

You could do a simple test to check if above two code snippets get the same result (Please check). One advantage of writing like this will make our code parallel over all the architechtures easily, we will see this soon.

Yes, this is due to the scheme in use. No need to worry about this now.

The behavior of the fluid is quite different from the dam_break_2d.py example however, and looks physically impossible especially along the left boundary. I put some pictures attached. Did I make some mistake, or is it related to the scheme used? (Vmagnitude is plotted).

Thanks!

We are getting closer in writing a simulation using equations. In the last assignment you wrote a script to compute the force on a set of particle where each particle influences each other particle. As an extension to this, suppose say we have two sets of particle arrays and we want to compute the force on each particle due all the particle in both the particle arrays.

A little context for you to start with.

1. Create two particle arrays with names as "bodies1" (with mass to be equal to all the particles) and "bodies2" (same conditions as of bodies but the position would be different).

2. If you know numerical integration of Newton's second law, try to integrate the system to the new time step. Then again compute the forces (here you have to make you forces zero before computation of the force) and integrate the system again.

Try to post a picture with the particle positions before the simulation and after the simulation (i.e., after two timesteps).

I have the next assignment ready, which also a simple extension as this, as soon as you submit.

[Enrique del Castillo]

Hi,

Q1a)  array                     Q2a) array               Q3a) error, size of array was specified and is incompatible

Q1b)  array                     Q2b) array               Q3b) error, same.

Q4a) error

Q4b) all values of length 10 array changed to 10.0

Below is my code for the two particle arrays problem, and I attached the figures as pictures. The particles seem to congregate around x=0, which seems a bit unusual. 

from pysph.base.particle_array import ParticleArray

from matplotlib import pyplot as plt

bodies1 = ParticleArray(name='fluid', x=[1,2,3,4,5,6,7,8,9,0], y=[.3,5,1,8,5,3,1,3,3,5],m=[0.1,2,1,3,4,5,6,7,2,4],k=[0.1,.3,.6,.35,.43,.5,.6,.7,.21,.9])

bodies2 = ParticleArray(name='fluid', x=[.5,2,4,3,.3,.56,2,7,.9,3], y=[3,7,2,.4,3,3,.1,8,9,4.5],m=[.3,.5,.2,1,5,4,.44,.32,.11,.9],k=[0.14,.22,.5,.25,.88,.94,.32,.17,.421,.24])

bodies1.add_property('tfx',type='double')

bodies1.add_property('tfy',type='double')

bodies1.add_property('vx',type='double')

bodies1.add_property('vy',type='double')

bodies1.add_property('ax',type='double')

bodies1.add_property('ay',type='double')

bodies2.add_property('tfx',type='double')

bodies2.add_property('tfy',type='double')

bodies2.add_property('vx',type='double')

bodies2.add_property('vy',type='double')

bodies2.add_property('ax',type='double')

bodies2.add_property('ay',type='double')

dt=.005

final_step = 2

plt.scatter(bodies1.x, bodies1.y, s=10, c='r', marker="o", label='bodies1')

plt.scatter(bodies2.x, bodies2.y, s=10, c='b', marker="o", label='bodies2')

plt.xlabel('x')

plt.ylabel('y')

plt.legend(loc='upper right')

plt.title('Original Particle Positions')

plt.show() 

for t in range(final_step):

    for i in range(len(bodies1.m)):

        bodies1.tfx[i] = 0

        bodies1.tfy[i] = 0

        bodies2.tfx[i] = 0

        bodies2.tfy[i] = 0

    #start with bodies1 particles            

    for i in range(len(bodies1.m)):

        # for a bodies1 particle consider interactions caused by all bodies1 particles

        for j in range(len(bodies1.m)): 

            if j!=i:

                k_i = bodies1.k[i]

                m_i = bodies1.m[i]

                m_j = bodies1.m[j]

                r_ij = ((bodies1.x[j]-bodies1.x[i])**2)+((bodies1.y[j]-bodies1.y[i])**2)**.5

                r_ij_2 = r_ij**2 

                F_ij = k_i*m_i*m_j/(r_ij_2)

                fx_ij = F_ij*((bodies1.x[j]-bodies1.x[i])/r_ij)

                fy_ij = F_ij*((bodies1.y[j]-bodies1.y[i])/r_ij)

                bodies1.tfx[i]+= fx_ij

                bodies1.tfy[i]+= fy_ij

        # for a bodies1 particle consider interactions caused by all bodies2 particles       

        for j in range(len(bodies2.m)):

            k_i = bodies1.k[i]

            m_i = bodies1.m[i]

            m_j = bodies2.m[j]

            r_ij = ((bodies2.x[j]-bodies1.x[i])**2)+((bodies2.y[j]-bodies1.y[i])**2)**.5

            r_ij_2 = r_ij**2 

            F_ij = k_i*m_i*m_j/(r_ij_2)

            fx_ij = F_ij*((bodies2.x[j]-bodies1.x[i])/r_ij)

            fy_ij = F_ij*((bodies2.y[j]-bodies1.y[i])/r_ij)

            bodies1.tfx[i]+= fx_ij

            bodies1.tfy[i]+= fy_ij   

    #now look at bodies2 particles       

    for i in range(len(bodies2.m)):

        # for a bodies2 particle consider interactions caused by all bodies2 particles

        for j in range(len(bodies2.m)):

            if j!=i:

                k_i = bodies2.k[i]

                m_i = bodies2.m[i]

                m_j = bodies2.m[j]

                r_ij = ((bodies2.x[j]-bodies2.x[i])**2)+((bodies2.y[j]-bodies2.y[i])**2)**.5

                r_ij_2 = r_ij**2 

                F_ij = k_i*m_i*m_j/(r_ij_2)

                fx_ij = F_ij*((bodies2.x[j]-bodies2.x[i])/r_ij)

                fy_ij = F_ij*((bodies2.y[j]-bodies2.y[i])/r_ij)

                bodies2.tfx[i]+= fx_ij

                bodies2.tfy[i]+= fy_ij

        # for a bodies2 particle consider interactions caused by all bodies1 particles       

        for j in range(len(bodies1.m)):

            k_i = bodies2.k[i]

            m_i = bodies2.m[i]

            m_j = bodies1.m[j]

            r_ij = ((bodies1.x[j]-bodies2.x[i])**2)+((bodies1.y[j]-bodies2.y[i])**2)**.5

            r_ij_2 = r_ij**2 

            F_ij = k_i*m_i*m_j/(r_ij_2)

            fx_ij = F_ij*((bodies1.x[j]-bodies2.x[i])/r_ij)

            fy_ij = F_ij*((bodies1.y[j]-bodies2.y[i])/r_ij)

            bodies2.tfx[i]+= fx_ij

            bodies2.tfy[i]+= fy_ij 

            

    # Update accelerations, velocities, and positions for bodies1 and bodies2 particles.               

    for i in range(len(bodies1.m)):

        bodies1.ax[i] = bodies1.tfx[i]/bodies1.m[i]

        bodies1.ay[i] = bodies1.tfy[i]/bodies1.m[i]

        bodies1.vx[i] = bodies1.vx[i] + dt*bodies1.ax[i]

        bodies1.vy[i] = bodies1.vy[i] + dt*bodies1.ay[i]

        bodies1.x[i] = bodies1.x[i] + dt*bodies1.vx[i]

        bodies1.y[i] = bodies1.y[i] + dt*bodies1.vy[i]

        bodies2.ax[i] = bodies2.tfx[i]/bodies2.m[i]

        bodies2.ay[i] = bodies2.tfy[i]/bodies2.m[i]

        bodies2.vx[i] = bodies2.vx[i] + dt*bodies2.ax[i]

        bodies2.vy[i] = bodies2.vy[i] + dt*bodies2.ay[i]

        bodies2.x[i] = bodies2.x[i] + dt*bodies2.vx[i]

        bodies2.y[i] = bodies2.y[i] + dt*bodies2.vy[i]

    plt.scatter(bodies1.x, bodies1.y, s=10, c='r', marker="o", label='bodies1')

    plt.scatter(bodies2.x, bodies2.y, s=10, c='b', marker="o", label='bodies2')

    plt.xlabel('x')

    plt.ylabel('y')

    plt.legend(loc='upper right')

    plt.title('Particle Positions after timestep: T=%i' %(t+1))

    plt.show() 

Thanks!

[A Dinesh]

On Tue, Feb 11, 2020 at 6:55 AM Enrique del Castillo <em...@alumni.princeton.edu> wrote:

Hi,

Q1a)  array                     Q2a) array               Q3a) error, size of array was specified and is incompatible

Q1b)  array                     Q2b) array               Q3b) error, same.

Q4a) error

Q4b) all values of length 10 array changed to 10.0

Below is my code for the two particle arrays problem, and I attached the figures as pictures. The particles seem to congregate around x=0, which seems a bit unusual. 

 

The force law I provided is just a prototype and physics is not to be concerned.

from pysph.base.particle_array import ParticleArray

from matplotlib import pyplot as plt

bodies1 = ParticleArray(name='fluid', x=[1,2,3,4,5,6,7,8,9,0], y=[.3,5,1,8,5,3,1,3,3,5],m=[0.1,2,1,3,4,5,6,7,2,4],k=[0.1,.3,.6,.35,.43,.5,.6,.7,.21,.9])

bodies2 = ParticleArray(name='fluid', x=[.5,2,4,3,.3,.56,2,7,.9,3], y=[3,7,2,.4,3,3,.1,8,9,4.5],m=[.3,.5,.2,1,5,4,.44,.32,.11,.9],k=[0.14,.22,.5,.25,.88,.94,.32,.17,.421,.24])

 

The name of the particle arrays should not be same.

bodies1.add_property('tfx',type='double')

bodies1.add_property('tfy',type='double')

bodies1.add_property('vx',type='double')

bodies1.add_property('vy',type='double')

bodies1.add_property('ax',type='double')

bodies1.add_property('ay',type='double')

bodies2.add_property('tfx',type='double')

bodies2.add_property('tfy',type='double')

bodies2.add_property('vx',type='double')

bodies2.add_property('vy',type='double')

bodies2.add_property('ax',type='double')

bodies2.add_property('ay',type='double')

You do the same with a function ,

```

    properties = [
        'uhat', 'vhat', 'what', 'rho0', 'rho_div', 'p0', 'auhat', 'avhat',
        'awhat', 'arho', 'arho0', 'cs', 'rho_tmp', 'rho_div'
    ]

    for i in properties:
        pa.add_property(i)

```

# Part 1

In the next part we will see how to use `PySPH` to execute the same equations on our particle arrays. Before that I want to introduce a new force law (a simple one).

```

fij = k_i * m_j * rij        # if rij < 1.0

fij = 0                          # if rij > 1.0

```

Where rij is the distance between the particles. The problem description is as follows, create particles  on a grid of length 1.0 and height 1.0,  with spacing 0.2.

Take the density to be 1000. Set mass accordingly. Create two such particle arrays one top on another with 0.1 gap between them. make sure particles don't fall on top)

In the previous two cases we were looping over all the particles to find the force (O(n^2)). Our new force law is limited and doesn't influence every other particle, we can use nearest neighbour

law and eliminate all the particles which are at certain distance, and only consider the particles which might exert a force. Based on this logic lets proceed to pysph.

In the previous assignment, we have two particle arrays, `bodies1` and `bodies2`, and every particle influences each other.

In short notation (I am writing psuedo code, for one dimension) the new force law can be written as

```

# to compute force on bodies1

# this is destination (implies, we want to find the force on bodies 1)

d_x = bodies1.x

d_ax = bodies1.ax

# Then what are the sources of bodies1?

# They are both, bodies 1 and bodies 2

s_x = bodies1.x

# now the loop goes as
for d_idx in range(len(d_x)):
    # where d_idx is the id of my destination particles (a notation used in PySPH)
    for s_idx in range(len(s_x)):
        # where s_idx is the id of my destination particles (a notation used in PySPH)
        rij = abs(s_x[s_idx])

        if rij > 1.0:
            d_ax[d_idx] += d_k[d_idx] * s_m[s_idx] * rij

# --------------------------

# repeat the same for bodies1 with sources as bodies2 ( Assignment )

# -------------------------

# Compute the force on bodes 2 (now bodies2 will be the destination, and bodies1, bodies 2 will the sources)

# try to integrate them too and plot the particle positions. (Please use K value to be same for all particles of a given particle array).

# This is simple and just by changing i to d_idx and j to s_idx should complete the assignment

```

Part 2

Even though the force law influence till some extent, we are not taking advantage of it and looping over all the particles (our algorithm is still O(n^2)).

As I discussed before one way to make our algorithms faster is to use nearest neighbours concept. The code,

```

for d_idx in range(len(d_x)):
    # get all the source particle in the influence range of d_idx
    neighbours = NNPS.get_neighbours(d_idx, source_name="bodies2")
    # Here NNPS is some black box, which will give us all the indices of the
    # source which will possibly influence d_idx
    # the function signature is not exact (is just for illustration)
    # neighbours is the list of s_idx (indices of sources possibly influencing d_idx)

    for s_idx in neighbours:
        # where s_idx is the id of my destination particles (a notation used in PySPH)
        rij = abs(s_x[s_idx])

        if rij > 1.0:
            d_ax[d_idx] += d_k[d_idx] * s_m[s_idx] * rij

```

Assignment:

Complete the above loop (psuedo) for other particle arrays (bodies1, bodies2).

Part 3

Before we implement this in terms of PySPH `Equation`, we will execute some PySPH equations as practice.

Create particles on grid using `get_particle_array_wcsph` function (link). Make sure you don't leave particles h, m, rho to be zero. 

Before we execute any equation ("here we will do summation density"), we will save the data of our particle array as follows

```

    # ----------------
    # DUMP THE DATA BEFORE RUNNING THE EQUATIONS
    # ----------------
    folder_name = "learn_running_equation_output"
    os.makedirs(folder_name, exist_ok=True)
    dump(os.path.join(filename, 'learn_running_equation_0'), [pa],
         dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=False)
    # ----------------
    # DUMP THE DATA BEFORE RUNNING THE EQUATIONS ENDS
    # ----------------

```

Where `pa` is your particle array. This will create a folder in the current directory. (Assignment)  To visualize this file do `pysph view learn_running_equation_output`.

Then execute the `SummationDensity` equation on our particle array as follows

```

    from pysph.base.kernels import QuinticSpline

    from pysph.sph.equation import Group
    from pysph.tools.sph_evaluator import SPHEvaluator

    from pysph.solver.output import dump

    from pysph.sph.basic_equations import SummationDensity

    # ---------------------
    # kernel, we will come to this later point (Once we know SPH, this will be clear)
    # ---------------------
    kernel = QuinticSpline(dim=dim)
    # ---------------------
    # kernel ends
    # ---------------------

    # ---------------------
    # Create equations
    # ---------------------
    pa_name = pa.name
    # create all the particles
    eqns = [
        Group(
            equations=[
                SummationDensity(dest=pa_name, sources=[pa_name]),
            ], )
    ]
    sph_eval = evaluate([pa], eqns, kernel, dim)
    sph_eval.evaluate()
    # ---------------------

    # ----------------
    # DUMP THE AFTER RUNNING THE EQUATIONS
    # ----------------

    dump(os.path.join(folder_name, 'learn_running_equation_1'), [pa],
         dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=True)
    # ----------------
    # DUMP THE DATA AFTER RUNNING THE EQUATIONS ENDS
    # ----------------

```

Now again do . (Assignment)   `pysph view learn_running_equation_output`. No you will find two files in the viewer (or folder).

Now check the density.

In the next part we will write our new force law in terms of PySPH equation format and we execute on our particle array and see the result

--
You received this message because you are subscribed to the Google Groups "pysph-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/7763cceb-6c00-4122-b054-5c616f092dd6%40googlegroups.com.

[Prabhu Ramachandran]

On Tuesday 11 February 2020 10:53 AM, A Dinesh wrote:

Where `pa` is your particle array. This will create a folder in the current directory. (Assignment)  To visualize this file do `pysph view learn_running_equation_output`.

Then execute the `SummationDensity` equation on our particle array as follows

    sph_eval = evaluate([pa], eqns, kernel, dim)

Should be

sph_eval = SPHEvaluator([pa], eqns, kernel, dim)

Prabhu

[Enrique del Castillo]

Hi, 

I've had some trouble with parts of the assignment.  For part 1 this is my code:

from __future__ import print_function

from pysph.base.particle_array import ParticleArray

import numpy as np

from matplotlib import pyplot as plt

dx = .2

dy = .2

x, y = np.mgrid[0:1.2:dx, 0:1.2:dy]

x2, y2 = np.mgrid[.1:1.3:dx, .1:1.3:dy]

bodies1 = ParticleArray(name='bodies1', x=x, y=y,m=1.0,k=.5)

bodies2 = ParticleArray(name='bodies2', x=x2, y=y2,m=1.0,k=.2)

plt.scatter(bodies1.x, bodies1.y, s=10, c='r', marker="o", label='bodies1')

plt.scatter(bodies2.x, bodies2.y, s=10, c='b', marker="o", label='bodies2')

plt.xlabel('x')

plt.ylabel('y')

plt.legend(loc='upper right')

plt.title('Original Particle Positions')

plt.show()

properties = ['tfx', 'tfy', 'vx', 'vy', 'ax', 'ay']

for i in properties:

    bodies1.add_property(i,type='double')

    bodies2.add_property(i,type='double')

    

dt=.005

final_step = 2

for t in range(final_step):

    

    for i in range(len(bodies1.m)):

        bodies1.tfx[i] = 0

        bodies1.tfy[i] = 0

        bodies2.tfx[i] = 0

        bodies2.tfy[i] = 0

    #compute force on bodies1

    #destination

    d_x = bodies1.x

    d_y = bodies1.y

    d_k = bodies1.k

    d_tfx = bodies1.tfx

    d_tfy = bodies1.tfy

    d_ax = bodies1.ax

    d_ay = bodies1.ay

    d_vx = bodies1.vx

    d_vy = bodies1.vy

    d_m = bodies1.m

    #source1

    s_x = bodies1.x

    s_y = bodies1.y

    s_m = bodies1.m

    for d_idx in range(len(d_x)):

        for s_idx in range(len(s_x)):

            r_ij = ((s_x[s_idx]-d_x[d_idx])**2)+((s_y[s_idx]-d_y[d_idx])**2)**.5

            #print(r_ij)

            if r_ij>1.0:

                F_ij = d_k[d_idx]*s_m[s_idx]*r_ij

            else:

                F_ij = 0.0

            fx_ij = F_ij*((s_x[s_idx]-d_x[d_idx])/r_ij)

            fy_ij = F_ij*((s_y[s_idx]-d_y[d_idx])/r_ij)

            d_tfx[d_idx]+= fx_ij

            d_tfy[d_idx]+= fy_ij

        #change source to bodies2

        s_x = bodies2.x

        s_y = bodies2.y

        s_m = bodies2.m

        for s_idx in range(len(s_x)):

            r_ij = ((s_x[s_idx]-d_x[d_idx])**2)+((s_y[s_idx]-d_y[d_idx])**2)**.5

            if r_ij>1.0:

                F_ij = d_k[d_idx]*s_m[s_idx]*r_ij

            else:

                F_ij = 0.0

            fx_ij = F_ij*((s_x[s_idx]-d_x[d_idx])/r_ij)

            fy_ij = F_ij*((s_y[s_idx]-d_y[d_idx])/r_ij)

            d_tfx[d_idx]+= fx_ij

            d_tfy[d_idx]+= fy_ij

        

        d_ax[d_idx] = d_tfx[d_idx]/d_m[d_idx]

        d_ay[d_idx] = d_tfy[d_idx]/d_m[d_idx]

        d_vx[d_idx] = d_vx[d_idx] + dt*d_ax[d_idx]

        d_vy[d_idx] = d_vy[d_idx] + dt*d_ay[d_idx]

        bodies1.tfx = d_tfx  

        bodies1.tfy = d_tfy 

        bodies1.ax = d_ax 

        bodies1.ay = d_ay  

        bodies1.vx = d_vx

        bodies1.vy = d_vy

        

    #compute force on bodies2

    #destination 

    d_x = bodies2.x

    d_y = bodies2.y

    d_k = bodies2.k

    d_tfx = bodies2.tfx

    d_tfy = bodies2.tfy

    d_ax = bodies2.ax

    d_ay = bodies2.ay

    d_vx = bodies2.vx

    d_vy = bodies2.vy

    d_m = bodies2.m

    #source1

    s_x = bodies2.x

    s_y = bodies2.y

    s_m = bodies2.m

    for d_idx in range(len(d_x)):

        for s_idx in range(len(s_x)):

            r_ij = ((s_x[s_idx]-d_x[d_idx])**2)+((s_y[s_idx]-d_y[d_idx])**2)**.5

            if r_ij>1.0:

                F_ij = d_k[d_idx]*s_m[s_idx]*r_ij

            else:

                F_ij = 0.0

            fx_ij = F_ij*((s_x[s_idx]-d_x[d_idx])/r_ij)

            fy_ij = F_ij*((s_y[s_idx]-d_y[d_idx])/r_ij)

            d_tfx[d_idx]+= fx_ij

            d_tfy[d_idx]+= fy_ij

        #change source to bodies2

        s_x = bodies1.x

        s_y = bodies1.y

        s_m = bodies1.m

        for s_idx in range(len(s_x)):

            r_ij = ((s_x[s_idx]-d_x[d_idx])**2)+((s_y[s_idx]-d_y[d_idx])**2)**.5

            if r_ij>1.0:

                F_ij = d_k[d_idx]*s_m[s_idx]*r_ij

            else:

                F_ij = 0.0

            fx_ij = F_ij*((s_x[s_idx]-d_x[d_idx])/r_ij)

            fy_ij = F_ij*((s_y[s_idx]-d_y[d_idx])/r_ij)

            d_tfx[d_idx]+= fx_ij

            d_tfy[d_idx]+= fy_ij

        d_ax[d_idx] = d_tfx[d_idx]/d_m[d_idx]

        d_ay[d_idx] = d_tfy[d_idx]/d_m[d_idx]

        d_vx[d_idx] = d_vx[d_idx] + dt*d_ax[d_idx]

        d_vy[d_idx] = d_vy[d_idx] + dt*d_ay[d_idx]

        bodies2.tfx = d_tfx  

        bodies2.tfy = d_tfy 

        bodies2.ax = d_ax 

        bodies2.ay = d_ay  

        bodies2.vx = d_vx

        bodies2.vy = d_vy

        #print(d_vx)

    # Update accelerations, velocities, and positions for bodies1 and bodies2 particles.               

    for i in range(len(s_x)):

        bodies1.x[i] = bodies1.x[i] + dt*bodies1.vx[i]

        bodies1.y[i] = bodies1.y[i] + dt*bodies1.vy[i]

        bodies2.x[i] = bodies2.x[i] + dt*bodies2.vx[i]

        bodies2.y[i] = bodies2.y[i] + dt*bodies2.vy[i]

    plt.scatter(bodies1.x, bodies1.y, s=10, c='r', marker="o", label='bodies1')

    plt.scatter(bodies2.x, bodies2.y, s=10, c='b', marker="o", label='bodies2')

    plt.xlabel('x')

    plt.ylabel('y')

    plt.legend(loc='upper right')

    plt.title('Particle Positions after timestep: T=%i' %(t+1))

    plt.show()        

I receive an error RuntimeWarning: invalid value encountered in double_scalars, and many of my d_tfx and d_tfy values become NaNs. 

For part 2) this is my loop: similar for the other cases changing the destination and sources.

    #example where destination is bodies1 and source is bodies2

    d_x = bodies1.x

    d_y = bodies1.y

    d_k = bodies1.k

    d_tfx = bodies1.tfx

    d_tfy = bodies1.tfy

    s_x = bodies1.x

    s_y = bodies1.y

    s_m = bodies1.m

    for d_idx in range(len(d_x)):

        neighbours = NNPS.get_neighbours(d_idx, source_name="bodies2")

        for s_idx in neighbors:

            r_ij = ((s_x[s_idx]-d_x[d_idx])**2)+((s_y[s_idx]-d_y[d_idx])**2)**.5

            if r_ij>1.0:

                F_ij = d_k[d_idx]*s_m[s_idx]*r_ij

            else:

                F_ij = 0.0

            fx_ij = F_ij*((s_x[s_idx]-d_x[d_idx])/r_ij)

            fy_ij = F_ij*((s_y[s_idx]-d_y[d_idx])/r_ij)

            d_tfx[d_idx]+= fx_ij

            d_tfy[d_idx]+= fy_ij

Now for Part 3) I can complete the first script, and visualize the particles (attached picture), but I cannot complete the next part. There seems to be an issue with the line

sph_eval = SPHEvaluator([pa], eqns, kernel, dim). I first needed to set dim to an integer, so I chose 2. But then I receive the error message AttributeError: 'int' object has no attribute '__dict__'

Thanks

[Enrique del Castillo]

forgot the picture

[A Dinesh]

Can you post the code for the third assignment?

--

You received this message because you are subscribed to the Google Groups "pysph-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/3d34e7bf-2ba7-47c7-adcb-60e7aae9d2bd%40googlegroups.com.

[Enrique del Castillo]

Sure this is what I had for the 3rd part. I was just using what you gave in the instructions and added u and v like in the elliptical drop example.

import os

import pickle

from pysph.base.utils import get_particle_array_wcsph

dx = .2

dy = .2

x, y = np.mgrid[0:1.2:dx, 0:1.2:dy]

m=3

h=.1

rho=10

u = -100*x

v = 100*y

name = 'test'

pa = get_particle_array_wcsph(x=x, y=y, m=m, rho=rho, h=h,u=u, v=v, name=name)

# ----------------

# DUMP THE DATA BEFORE RUNNING THE EQUATIONS

# ----------------

folder_name = "learn_running_equation_output3"

os.makedirs(folder_name)

dump(os.path.join(folder_name, 'learn_running_equation_0'), [pa],dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=False)

# ----------------

# DUMP THE DATA BEFORE RUNNING THE EQUATIONS ENDS

# ----------------

from pysph.base.kernels import QuinticSpline

    from pysph.sph.equation import Group

    from pysph.tools.sph_evaluator import SPHEvaluator

    from pysph.solver.output import dump

    from pysph.sph.basic_equations import SummationDensity

    dim=2

    # ---------------------

    # kernel, we will come to this later point (Once we know SPH, this will be clear)

    # ---------------------

    kernel = QuinticSpline(dim=2)

    # ---------------------

    # kernel ends

    # ---------------------

    # ---------------------

    # Create equations

    # ---------------------

    pa_name = pa.name

    # create all the particles

    eqns = [

        Group(

            equations=[

                SummationDensity(dest=pa_name, sources=[pa_name]),

            ], )

    ]

    sph_eval = SPHEvaluator([pa], eqns, kernel, dim)

    sph_eval.evaluate()

    # ---------------------

    # ----------------

    # DUMP THE AFTER RUNNING THE EQUATIONS

    # ----------------

    dump(os.path.join(folder_name, 'learn_running_equation_1'), [pa],

         dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=True)

    # ----------------

    # DUMP THE DATA AFTER RUNNING THE EQUATIONS ENDS

    # ----------------

To unsubscribe from this group and stop receiving emails from it, send an email to pysph...@googlegroups.com.

[A Dinesh]

I see the problem. It should be

```

sph_eval = SPHEvaluator([pa], equations=eqns, kernel=kernel, dim=dim)

```

Use high resolution (make dx small). Try looking at the density after running the equations.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/5c0e9a85-8135-4a13-9c98-e965e85c60ba%40googlegroups.com.

[Enrique del Castillo]

When I try running this corrected version in the terminal I receive the error: 

Traceback (most recent call last):

  File "density.py", line 25, in <module>

    dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=False)

TypeError: dump() takes at most 4 arguments (5 given)

while if I try to run it in a jupyter notebook (where I only have python2.7)

I receive the warning:

warning: include path for stdlibc++ headers not found; pass '-stdlib=libc++' on the command line to use the libc++ standard library instead [-Wstdlibcxx-not-found]

and the error:

error: 'ios' file not found
#include "ios"

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/5c0e9a85-8135-4a13-9c98-e965e85c60ba%40googlegroups.com.

[A Dinesh]

Can you try this?

```

import os
import pickle
from pysph.base.utils import get_particle_array_wcsph

from pysph.base.kernels import QuinticSpline
from pysph.sph.equation import Group
from pysph.tools.sph_evaluator import SPHEvaluator
from pysph.solver.output import dump
from pysph.sph.basic_equations import SummationDensity

import numpy as np

dx = .2
dy = .2
x, y = np.mgrid[0:1.2:dx, 0:1.2:dy]
m=3
h=.1
rho=10
u = -100*x
v = 100*y
name = 'test'
pa = get_particle_array_wcsph(x=x, y=y, m=m, rho=rho, h=h,u=u, v=v, name=name)

# ----------------
# DUMP THE DATA BEFORE RUNNING THE EQUATIONS
# ----------------
folder_name = "learn_running_equation_output3"

os.makedirs(folder_name, exist_ok=True)

dump(os.path.join(folder_name, 'learn_running_equation_0'), [pa],dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=False)
# ----------------
# DUMP THE DATA BEFORE RUNNING THE EQUATIONS ENDS
# ----------------

dim=2
# ---------------------
# kernel, we will come to this later point (Once we know SPH, this will be clear)
# ---------------------
kernel = QuinticSpline(dim=2)
# ---------------------
# kernel ends
# ---------------------

# ---------------------
# Create equations
# ---------------------
pa_name = pa.name
# create all the particles
eqns = [
    Group(
        equations=[
            SummationDensity(dest=pa_name, sources=[pa_name]),
        ], )
]

sph_eval = SPHEvaluator([pa], equations=eqns, kernel=kernel, dim=dim)

sph_eval.evaluate()
# ---------------------

# ----------------
# DUMP THE AFTER RUNNING THE EQUATIONS
# ----------------
dump(os.path.join(folder_name, 'learn_running_equation_1'), [pa],
        dict(t=0, dt=0.1, count=0), detailed_output=True, only_real=True)
# ----------------
# DUMP THE DATA AFTER RUNNING THE EQUATIONS ENDS
# ----------------

```

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/16466983-127a-44e3-ab70-708da78ac138%40googlegroups.com.

[Enrique del Castillo]

I still get the same error unfortunately. Note I had to change the makedirs() to not take the exist_ok argument as I was also getting this error before.

makedirs() got an unexpected keyword argument 'exist_ok'

Is it possible to continue to the next step perhaps?

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/16466983-127a-44e3-ab70-708da78ac138%40googlegroups.com.

[A Dinesh]

We can go ahead, but make sure you understand this part. Try to run these in `python3` not `python2`.

As part of the next step, we need to read some documentation on how `PySPH` implements equations. In hand on tasks we were repeating ourselves, to find the force on a particle

due to several other particle arrays. By using the `Equation` class we were able to make it simple and reuse. Now we will look at using Groups and understand how `PySPH` eliminates

recomputation of several predicates before finding the force.

https://pysph.readthedocs.io/en/latest/design/equations.html

After this part we will implement basic rigid bodies colliding among them selves.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/e43ed226-97bf-49ae-a5a2-27451ba7a7f4%40googlegroups.com.

[A Dinesh]

Specifically read this section,

https://pysph.readthedocs.io/en/latest/design/equations.html#understanding-groups-a-bit-more

[Enrique del Castillo]

Ok I read through the equations documentation and looked at the examples in the SPH equations page. It makes sense to me.

On Tuesday, February 18, 2020 at 9:04:05 AM UTC-8, A Dinesh wrote:

Specifically read this section,

https://pysph.readthedocs.io/en/latest/design/equations.html#understanding-groups-a-bit-more

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/e43ed226-97bf-49ae-a5a2-27451ba7a7f4%40googlegroups.com.

[A Dinesh]

We will look at simulating rigid bodies in PySPH with the help of bouncing cube example. You can find the full code at https://github.com/pypr/pysph/blob/master/pysph/examples/rigid_body/bouncing_cube.py.

To explain,

"""A cube bouncing inside a box. (5 seconds)

This is used to test the rigid body equations.

"""

import numpy as np

from pysph.base.kernels import CubicSpline

from pysph.base.utils import get_particle_array_rigid_body

from pysph.sph.equation import Group

from pysph.sph.integrator import EPECIntegrator

from pysph.solver.application import Application

from pysph.solver.solver import Solver

Import the requisites for rigid body simulation.

from pysph.sph.rigid_body import (BodyForce, RigidBodyCollision,

RigidBodyMoments, RigidBodyMotion,

RK2StepRigidBody)

dim = 3

dt = 5e-3

tf = 5.0

gz = -9.81

hdx = 1.0

dx = dy = 0.02

rho0 = 10.0

class BouncingCube(Application):

def create_particles(self):

nx, ny, nz = 10, 10, 10

dx = 1.0 / (nx - 1)

x, y, z = np.mgrid[0:1:nx * 1j, 0:1:ny * 1j, 0:1:nz * 1j]

x = x.flat

y = y.flat

z = (z - 1).flat

m = np.ones_like(x) * dx * dx * rho0

h = np.ones_like(x) * hdx * dx

# radius of each sphere constituting in cube

rad_s = np.ones_like(x) * dx

Create the rigid body

body = get_particle_array_rigid_body(name='body', x=x, y=y, z=z, h=h,

m=m, rad_s=rad_s)

Set its center of mass velocity

body.vc[0] = -5.0

body.vc[2] = -5.0

# Create the tank.

nx, ny, nz = 40, 40, 40

dx = 1.0 / (nx - 1)

xmin, xmax, ymin, ymax, zmin, zmax = -2, 2, -2, 2, -2, 2

x, y, z = np.mgrid[xmin:xmax:nx * 1j, ymin:ymax:ny * 1j, zmin:zmax:nz *

1j]

interior = ((x < 1.8) & (x > -1.8)) & ((y < 1.8) & (y > -1.8)) & (

(z > -1.8) & (z <= 2))

tank = np.logical_not(interior)

x = x[tank].flat

y = y[tank].flat

z = z[tank].flat

m = np.ones_like(x) * dx * dx * rho0

h = np.ones_like(x) * hdx * dx

# radius of each sphere constituting in cube

rad_s = np.ones_like(x) * dx

tank = get_particle_array_rigid_body(name='tank', x=x, y=y, z=z, h=h,

m=m, rad_s=rad_s)

tank.total_mass[0] = np.sum(m)

return [body, tank]

def create_solver(self):

kernel = CubicSpline(dim=dim)

Use EPEC integrator for rigid body update

integrator = EPECIntegrator(body=RK2StepRigidBody())

solver = Solver(kernel=kernel, dim=dim, integrator=integrator, dt=dt,

tf=tf, adaptive_timestep=False)

solver.set_print_freq(10)

return solver

def create_equations(self):

equations = [

Group(equations=[

Apply body force for the body

BodyForce(dest='body', sources=None, gz=gz),

Resolve the interaction between the rigid body and tank with RigidBodyCollision equation.

RigidBodyCollision(dest='body', sources=['tank'], kn=1e4, en=1)

]),

After finding the forces, find the moment at the center of mass due to the forces acting on the particles constitued by the body.

Group(equations=[RigidBodyMoments(dest='body', sources=None)]),

Find the velocity of the particles from center of mass linear and angualr velocity.

Group(equations=[RigidBodyMotion(dest='body', sources=None)]),

]

return equations

if __name__ == '__main__':

app = BouncingCube()

Run the application.

app.run()

This should get you started. As part of the assignment, try to creating a few rigid bodies (May be 2 or 3) with arbitrary shape and set their  center of mass velocities such that

they will collide. Modify the equations so that thier collision is resolved.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/466ba586-6e65-4277-b746-68a2184c7235%40googlegroups.com.

[Enrique del Castillo]

Ok, I understand this example. The only issue is when I try to run it, I receive the same error I was getting before, 

ERROR

running build_ext

fatal error: 'ios' file not found

#include "ios"

and the warnings:

warning: include path for stdlibc++ headers not found; pass '-stdlib=libc++' on the command line to use the libc++ standard library instead [-Wstdlibcxx-not-found]

#warning "Using deprecated NumPy API, disable it with " \

Might you have any idea what is causing this?

Thanks!

<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border

[Prabhu Ramachandran]

On 20/02/20 12:16 PM, Enrique del Castillo wrote:

Ok, I understand this example. The only issue is when I try to run it, I receive the same error I was getting before, 

ERROR

running build_ext

fatal error: 'ios' file not found

#include "ios"

and the warnings:

warning: include path for stdlibc++ headers not found; pass '-stdlib=libc++' on the command line to use the libc++ standard library instead [-Wstdlibcxx-not-found]

#warning "Using deprecated NumPy API, disable it with " \

Again, it is likely because it is not finding the right compiler.  On a jupyter notebook you will need to set the CXX env var (see here: https://stackoverflow.com/questions/37890898/how-to-set-env-variable-in-jupyter-notebook).  There are many ways of doing this and some may work and some may not.  Sorry but jupyter notebooks do make some of this a bit more complex.  Regardless, I am not sure why clang is not installed properly with the c++ libraries?  Have you installed xcode?

Prabhu

[Enrique del Castillo]

Ok I managed to solve the problem. I had to reset the CXX env variable again, and it runs. I had already done this strangely, but it worked. I was using the terminal this time around too.

Anyways, ready to continue!

[A Dinesh]

Have you tried to tweak the above code a few more rigid bodies? Please try to change the above example with 3 rigid bodies and a tank, and allow them to collide each other. Run for 1 sec, and post the picture. Good luck.

--
You received this message because you are subscribed to the Google Groups "pysph-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/1ad96895-cd7b-423a-9b4e-a29ed9c40e67%40googlegroups.com.

[Enrique del Castillo]

I tried doing this, but for some reason I cannot keep 2/3 of the cubes from penetrating the box. I created two additional cubes, and it is these additional ones (body1 and body2) which are going through the boundary. I attached a picture. Also here is my code:

import numpy as np

from pysph.base.kernels import CubicSpline
from pysph.base.utils import get_particle_array_rigid_body
from pysph.sph.equation import Group

from pysph.sph.integrator import EPECIntegrator

from pysph.solver.application import Application
from pysph.solver.solver import Solver

from pysph.sph.rigid_body import (BodyForce, RigidBodyCollision,
                                  RigidBodyMoments, RigidBodyMotion,
                                  RK2StepRigidBody)

dim = 3

dt = 5e-3

tf = 1.0

gz = -9.81

hdx = 1.0
dx = dy = 0.02
rho0 = 10.0

class TWOcube(Application):

    def create_particles(self):
        nx, ny, nz = 10, 10, 10
        dx = 1.0 / (nx - 1)
        x, y, z = np.mgrid[0:1:nx * 1j, 0:1:ny * 1j, 0:1:nz * 1j]

        x1, y1, z1 = np.mgrid[.5:1.5:nx * 1j, .5:1.5:ny * 1j, .5:1.5:nz * 1j]
        x2, y2, z2 = np.mgrid[.5:1.5:nx * 1j, .5:1.5:ny * 1j, .5:1.5:nz * 1j]

        x = x.flat
        y = y.flat
        z = (z - 1).flat

        x1 = x1.flat
        y1 = y1.flat
        z1 = (z1).flat
        x2 = x2.flat
        y2 = y2.flat
        z2 = (z2 - 1.2).flat

        m = np.ones_like(x) * dx * dx * rho0
        h = np.ones_like(x) * hdx * dx
        # radius of each sphere constituting in cube
        rad_s = np.ones_like(x) * dx

        body = get_particle_array_rigid_body(name='body', x=x, y=y, z=z, h=h,
                                             m=m, rad_s=rad_s)

        body1 = get_particle_array_rigid_body(name='body1', x=x1, y=y1, z=z1, h=h,
                                              m=m, rad_s=rad_s)
        body2 = get_particle_array_rigid_body(name='body2', x=x2, y=y2, z=z2, h=h,
                                              m=m, rad_s=rad_s)

        body.vc[0] = -5.0
        body.vc[2] = -5.0

        body1.vc[0] = -5.0
        body1.vc[2] = -5.0
        body2.vc[0] = -5.0
        body2.vc[2] = -5.0

        # Create the tank.
        nx, ny, nz = 40, 40, 40
        dx = 1.0 / (nx - 1)
        xmin, xmax, ymin, ymax, zmin, zmax = -2, 2, -2, 2, -2, 2
        x, y, z = np.mgrid[xmin:xmax:nx * 1j, ymin:ymax:ny * 1j, zmin:zmax:nz *
                                                                           1j]
        interior = ((x < 1.8) & (x > -1.8)) & ((y < 1.8) & (y > -1.8)) & (
                (z > -1.8) & (z <= 2))
        tank = np.logical_not(interior)
        x = x[tank].flat
        y = y[tank].flat
        z = z[tank].flat
        m = np.ones_like(x) * dx * dx * rho0
        h = np.ones_like(x) * hdx * dx

        # radius of each sphere constituting in cube
        rad_s = np.ones_like(x) * dx
        tank = get_particle_array_rigid_body(name='tank', x=x, y=y, z=z, h=h,
                                             m=m, rad_s=rad_s)
        tank.total_mass[0] = np.sum(m)

        return [body2, body, body1, tank]

    def create_solver(self):
        kernel = CubicSpline(dim=dim)

        integrator = EPECIntegrator(body=RK2StepRigidBody(), body1=RK2StepRigidBody(), body2=RK2StepRigidBody())

        solver = Solver(kernel=kernel, dim=dim, integrator=integrator, dt=dt,
                        tf=tf, adaptive_timestep=False)
        solver.set_print_freq(10)
        return solver

    def create_equations(self):
        equations = [
            Group(equations=[

                BodyForce(dest='body', sources=None, gz=gz),

                RigidBodyCollision(dest='body', sources=['tank', 'body1', 'body2'], kn=1e4, en=1)
            ]),
            Group(equations=[
                BodyForce(dest='body1', sources=None, gz=gz),
                RigidBodyCollision(dest='body1', sources=['tank', 'body', 'body2'], kn=1e4, en=1)
            ]),
            Group(equations=[
                BodyForce(dest='body2', sources=None, gz=gz),
                RigidBodyCollision(dest='body2', sources=['tank', 'body1', 'body'], kn=1e4, en=1)
            ]),

            Group(equations=[RigidBodyMoments(dest='body', sources=None)]),

            Group(equations=[RigidBodyMotion(dest='body', sources=None)]),

            Group(equations=[RigidBodyMoments(dest='body1', sources=None)]),
            Group(equations=[RigidBodyMotion(dest='body1', sources=None)]),
            Group(equations=[RigidBodyMoments(dest='body2', sources=None)]),
            Group(equations=[RigidBodyMotion(dest='body2', sources=None)]),

        ]
        return equations


if __name__ == '__main__':

    app = TWOcube()
    app.run()

Is it an issue with the create_equations method?

Thanks!

To unsubscribe from this group and stop receiving emails from it, send an email to pysph...@googlegroups.com.

[A Dinesh]

Code looks OK. Try increasing the stiffness (kn) so that particles won't penetrate. If you Increase the stiffness of spring then you should also decrease the time step.

Or try increasing the spacing between the particles of the body.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/86e2b801-914d-4007-ab27-43ee1762f646%40googlegroups.com.

[Enrique del Castillo]

Thanks, increasing the stiffness, decreasing the timestep a bit and reducing the number of particles in the body did the trick in the end. I attached a picture after 1 second of simulation.

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/86e2b801-914d-4007-ab27-43ee1762f646%40googlegroups.com.

[A Dinesh]

As a next step, we will simulate too many rigid bodies. Before proceeding to write such, let us assume a problem where we have 20 rigid bodies, each have their name as

```

for i in range(30):

    body_i.name = "i"

```

That loop is a little funny, but I guess you got the point. Now inorder to write the interaction between each body, we will have, 30 equations, and 30 equations for finding the rigid body moments, etc.

But as the number of bodies increase, say 1000, we can't write so many equations (we can write them, but looks ugly). Inorder to overcome this problem, we create all the rigid bodies in a single particle

array, and give them body_id property. And the interaction is trivial to write.

As an example look at this example. Mainly this point https://github.com/pypr/pysph/blob/master/pysph/examples/rigid_body/bouncing_cubes.py#L68.

Now adapt your previous example to this formulation. Then create 100 bodies and run the simulation.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/4251d389-7f9a-4969-af21-d096162e5cce%40googlegroups.com.

[Enrique del Castillo]

I was able to write a script similar to bouncing_cubes.py keeping all the bodies within 1 particle array, but for 100 bodies, the simulation takes a long time to run, and after a couple of hours quits with the error 

RuntimeError: ERROR: LinkedListNNPS requires too many cells (268588224).

Was it expected to take so long to run?

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/4251d389-7f9a-4969-af21-d096162e5cce%40googlegroups.com.

[A Dinesh]

You will understand the possible reason by looking at the viewer. I think the bodies are blowing up. Or may be you could just try with 10 bodies. And run for a while. This part of the assignment is to mainly understand concept of body_id.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/76e0ffa0-5dad-4050-8564-fc4a38a998fb%40googlegroups.com.

[Enrique del Castillo]

Got it. Ok Here it is with 10 bodies in the screenshot. Ready for the next part.

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/76e0ffa0-5dad-4050-8564-fc4a38a998fb%40googlegroups.com.

[A Dinesh]

Till now we are able to solve the rigid body collisions and rigid body dynamics in PySPH. Lets now see how we can simulate fluid flow in PySPH. Before going ahead and start writing the examples with equations, lets use a few schemes

and run a few problems where there is only fluid (Specially free surface). 

Go to pysph examples directory and take some examples which have fluid (dam_break_2d or hydrostatic_tank or dam_break_3d). Try running it with different schemes. While you run the same example with different scheme, you can output

it in different folder for different schemes by using the following command

Assignment1:

```

python dam_break_2d.py --scheme aha -d dam_break_2d_aha_output

```

If you don't give `-d dam_break_2d_aha_output` then it will output in `dam_break_2d_output`. While you again run it with different scheme then those output files will also clutter the same directory. So keep different folders for different schemes. 

Assignment 2:

What files have you found in the output directory. Just report back the file types.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/6ccce229-4f7b-435e-85b0-c9364b08f263%40googlegroups.com.

[Enrique del Castillo]

Ok I ran a few examples changing the scheme and the output folder.

And for part 2) these are the files one gets in the output folder for scheme aha. There are .info, .log, .hdf5, and a mayav_config.py file.

dam_break_2d.info dam_break_2d_16000.hdf5 dam_break_2d_3700.hdf5

dam_break_2d.log dam_break_2d_16100.hdf5 dam_break_2d_3800.hdf5

dam_break_2d_0.hdf5 dam_break_2d_16200.hdf5 dam_break_2d_3900.hdf5

dam_break_2d_100.hdf5 dam_break_2d_16300.hdf5 dam_break_2d_400.hdf5

dam_break_2d_1000.hdf5 dam_break_2d_16400.hdf5 dam_break_2d_4000.hdf5

dam_break_2d_10000.hdf5 dam_break_2d_16500.hdf5 dam_break_2d_4100.hdf5

dam_break_2d_10100.hdf5 dam_break_2d_16600.hdf5 dam_break_2d_4200.hdf5

dam_break_2d_10200.hdf5 dam_break_2d_16700.hdf5 dam_break_2d_4300.hdf5

dam_break_2d_10280.hdf5 dam_break_2d_16800.hdf5 dam_break_2d_4400.hdf5

dam_break_2d_10300.hdf5 dam_break_2d_16900.hdf5 dam_break_2d_4500.hdf5

dam_break_2d_10400.hdf5 dam_break_2d_1700.hdf5 dam_break_2d_4600.hdf5

dam_break_2d_10500.hdf5 dam_break_2d_17000.hdf5 dam_break_2d_4700.hdf5

dam_break_2d_10600.hdf5 dam_break_2d_17100.hdf5 dam_break_2d_4800.hdf5

dam_break_2d_10700.hdf5 dam_break_2d_17200.hdf5 dam_break_2d_4900.hdf5

dam_break_2d_10800.hdf5 dam_break_2d_17300.hdf5 dam_break_2d_500.hdf5

dam_break_2d_10900.hdf5 dam_break_2d_17400.hdf5 dam_break_2d_5000.hdf5

dam_break_2d_1100.hdf5 dam_break_2d_17500.hdf5 dam_break_2d_5100.hdf5

dam_break_2d_11000.hdf5 dam_break_2d_17600.hdf5 dam_break_2d_5140.hdf5

dam_break_2d_11100.hdf5 dam_break_2d_17700.hdf5 dam_break_2d_5200.hdf5

dam_break_2d_11200.hdf5 dam_break_2d_17800.hdf5 dam_break_2d_5300.hdf5

dam_break_2d_11300.hdf5 dam_break_2d_17900.hdf5 dam_break_2d_5400.hdf5

dam_break_2d_11400.hdf5 dam_break_2d_1800.hdf5 dam_break_2d_5500.hdf5

dam_break_2d_11500.hdf5 dam_break_2d_18000.hdf5 dam_break_2d_5600.hdf5

dam_break_2d_11600.hdf5 dam_break_2d_18100.hdf5 dam_break_2d_5700.hdf5

dam_break_2d_11700.hdf5 dam_break_2d_18200.hdf5 dam_break_2d_5800.hdf5

dam_break_2d_11800.hdf5 dam_break_2d_18300.hdf5 dam_break_2d_5900.hdf5

dam_break_2d_11900.hdf5 dam_break_2d_18400.hdf5 dam_break_2d_600.hdf5

dam_break_2d_1200.hdf5 dam_break_2d_18500.hdf5 dam_break_2d_6000.hdf5

dam_break_2d_12000.hdf5 dam_break_2d_18600.hdf5 dam_break_2d_6100.hdf5

dam_break_2d_12100.hdf5 dam_break_2d_18700.hdf5 dam_break_2d_6200.hdf5

dam_break_2d_12200.hdf5 dam_break_2d_18800.hdf5 dam_break_2d_6300.hdf5

dam_break_2d_12300.hdf5 dam_break_2d_18900.hdf5 dam_break_2d_6400.hdf5

dam_break_2d_12400.hdf5 dam_break_2d_1900.hdf5 dam_break_2d_6500.hdf5

dam_break_2d_12500.hdf5 dam_break_2d_19000.hdf5 dam_break_2d_6600.hdf5

dam_break_2d_12600.hdf5 dam_break_2d_19100.hdf5 dam_break_2d_6700.hdf5

dam_break_2d_12700.hdf5 dam_break_2d_19200.hdf5 dam_break_2d_6800.hdf5

dam_break_2d_12800.hdf5 dam_break_2d_19300.hdf5 dam_break_2d_6900.hdf5

dam_break_2d_12850.hdf5 dam_break_2d_19400.hdf5 dam_break_2d_700.hdf5

dam_break_2d_12900.hdf5 dam_break_2d_19500.hdf5 dam_break_2d_7000.hdf5

dam_break_2d_1300.hdf5 dam_break_2d_19600.hdf5 dam_break_2d_7100.hdf5

dam_break_2d_13000.hdf5 dam_break_2d_19700.hdf5 dam_break_2d_7200.hdf5

dam_break_2d_13100.hdf5 dam_break_2d_19800.hdf5 dam_break_2d_7300.hdf5

dam_break_2d_13200.hdf5 dam_break_2d_19900.hdf5 dam_break_2d_7400.hdf5

dam_break_2d_13300.hdf5 dam_break_2d_200.hdf5 dam_break_2d_7500.hdf5

dam_break_2d_13400.hdf5 dam_break_2d_2000.hdf5 dam_break_2d_7600.hdf5

dam_break_2d_13500.hdf5 dam_break_2d_20000.hdf5 dam_break_2d_7700.hdf5

dam_break_2d_13600.hdf5 dam_break_2d_20100.hdf5 dam_break_2d_7710.hdf5

dam_break_2d_13700.hdf5 dam_break_2d_20200.hdf5 dam_break_2d_7800.hdf5

dam_break_2d_13800.hdf5 dam_break_2d_20300.hdf5 dam_break_2d_7900.hdf5

dam_break_2d_13900.hdf5 dam_break_2d_20400.hdf5 dam_break_2d_800.hdf5

dam_break_2d_1400.hdf5 dam_break_2d_20500.hdf5 dam_break_2d_8000.hdf5

dam_break_2d_14000.hdf5 dam_break_2d_20600.hdf5 dam_break_2d_8100.hdf5

dam_break_2d_14100.hdf5 dam_break_2d_20700.hdf5 dam_break_2d_8200.hdf5

dam_break_2d_14200.hdf5 dam_break_2d_20800.hdf5 dam_break_2d_8300.hdf5

dam_break_2d_14300.hdf5 dam_break_2d_20900.hdf5 dam_break_2d_8400.hdf5

dam_break_2d_14400.hdf5 dam_break_2d_2100.hdf5 dam_break_2d_8500.hdf5

dam_break_2d_14500.hdf5 dam_break_2d_21000.hdf5 dam_break_2d_8600.hdf5

dam_break_2d_14600.hdf5 dam_break_2d_2200.hdf5 dam_break_2d_8700.hdf5

dam_break_2d_14700.hdf5 dam_break_2d_2300.hdf5 dam_break_2d_8800.hdf5

dam_break_2d_14800.hdf5 dam_break_2d_2400.hdf5 dam_break_2d_8900.hdf5

dam_break_2d_14900.hdf5 dam_break_2d_2500.hdf5 dam_break_2d_900.hdf5

dam_break_2d_1500.hdf5 dam_break_2d_2600.hdf5 dam_break_2d_9000.hdf5

dam_break_2d_15000.hdf5 dam_break_2d_2700.hdf5 dam_break_2d_9100.hdf5

dam_break_2d_15100.hdf5 dam_break_2d_2800.hdf5 dam_break_2d_9200.hdf5

dam_break_2d_15200.hdf5 dam_break_2d_2900.hdf5 dam_break_2d_9300.hdf5

dam_break_2d_15300.hdf5 dam_break_2d_300.hdf5 dam_break_2d_9400.hdf5

dam_break_2d_15400.hdf5 dam_break_2d_3000.hdf5 dam_break_2d_9500.hdf5

dam_break_2d_15500.hdf5 dam_break_2d_3100.hdf5 dam_break_2d_9600.hdf5

dam_break_2d_15600.hdf5 dam_break_2d_3200.hdf5 dam_break_2d_9700.hdf5

dam_break_2d_15700.hdf5 dam_break_2d_3300.hdf5 dam_break_2d_9800.hdf5

dam_break_2d_15800.hdf5 dam_break_2d_3400.hdf5 dam_break_2d_9900.hdf5

dam_break_2d_15900.hdf5 dam_break_2d_3500.hdf5 mayavi_config.py

dam_break_2d_1600.hdf5 dam_break_2d_3600.hdf5

To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/6ccce229-4f7b-435e-85b0-c9364b08f263%40googlegroups.com.

[A Dinesh]

Try also with a few other schemes. You can find all the supported schemes with `python filename.py -h | grep scheme`.

After running with like 2 to 3 schemes, check the log file in each folder. As an assignment, report back the equations used in each scheme for a specific example.

To unsubscribe from this group and stop receiving emails from it, send an email to pysph-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pysph-users/7ae080d4-3ac5-40e4-8036-20a9bc84872e%40googlegroups.com.