Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
258 changes: 125 additions & 133 deletions src/amuse/ext/protodisk.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import numpy
import warnings
from amuse.support.helpers import rename_fn_par
from amuse.ext.evrard_test import body_centered_grid_unit_cube
from amuse.ext.evrard_test import regular_grid_unit_cube
from amuse.ext.evrard_test import uniform_random_unit_cube
Expand All @@ -8,157 +9,148 @@

from amuse.datamodel import Particles
from amuse.datamodel import ParticlesWithUnitsConverted


def approximate_inverse_error_function(x):
a=8*(numpy.pi-3)/3*numpy.pi*(4-numpy.pi)
return numpy.sign(x)*numpy.sqrt(
numpy.sqrt((2/numpy.pi/a+numpy.log(1-x**2)/2)**2-numpy.log(1-x**2)/a)-(2/numpy.pi/a+numpy.log(1-x**2)/2)
)

class uniform_unit_cylinder(object):
def __init__(self,targetN, base_grid=None):
cube_cylinder_ratio=numpy.pi*0.5**2
self.targetN=targetN
self.estimatedN=targetN/cube_cylinder_ratio
a = 8 * (numpy.pi - 3) / 3 * numpy.pi * (4 - numpy.pi)

return numpy.sign(x) * numpy.sqrt(
numpy.sqrt(
(2 / numpy.pi / a + numpy.log(1 - x**2) / 2)**2 - numpy.log(1 - x**2) / a
) - (2 / numpy.pi / a + numpy.log(1 - x**2) / 2)
)

class uniform_unit_cylinder:
def __init__(self, targetN, base_grid=None):
cube_cylinder_ratio = numpy.pi*0.5**2
self.targetN = targetN
self.estimatedN = targetN / cube_cylinder_ratio

if base_grid is None:
self.base_grid=uniform_random_unit_cube
self.base_grid = uniform_random_unit_cube
else:
self.base_grid=base_grid
def cutout_cylinder(self,x,y,z):
r=x**2+y**2
selection=r < numpy.ones_like(r)
x=x.compress(selection)
y=y.compress(selection)
z=z.compress(selection)
return x,y,z
self.base_grid = base_grid

def cutout_cylinder(self, x, y, z):
r = x**2 + y**2
selection = r < numpy.ones_like(r)
x = x.compress(selection)
y = y.compress(selection)
z = z.compress(selection)
return x, y, z

def make_xyz(self):
if(self.base_grid==uniform_random_unit_cube):
estimatedN=self.estimatedN
x=[]
if (self.base_grid == uniform_random_unit_cube):
estimatedN = self.estimatedN
x = []
while len(x) < self.targetN:
estimadedN=estimatedN*1.1+1
x,y,z=self.cutout_cylinder(*(self.base_grid(estimatedN)).make_xyz())
return x[0:self.targetN],y[0:self.targetN],z[0:self.targetN]
estimadedN = estimatedN * 1.1 + 1
x, y, z = self.cutout_cylinder(*(self.base_grid(estimatedN)).make_xyz())

return x[0:self.targetN], y[0:self.targetN], z[0:self.targetN]
else:
return self.cutout_cylinder(*(self.base_grid(self.estimatedN)).make_xyz())


class ProtoPlanetaryDisk:

def __init__(
self, targetN, convert_nbody=None, discfraction=0.1,
densitypower=1., thermalpower=0.5, radius_min=1, radius_max=100,
gamma=1.,q_out=2.,base_grid=None, Rmin=None, Rmax=None,
densitypower=1., thermalpower=0.5, radius_min=None, radius_max=None,
gamma=1., q_out=2., base_grid=None, Rmin=None, Rmax=None,
):
if Rmin is not None:
warnings.warn(
"Rmin is deprecated, use radius_min instead",
category=FutureWarning,
)
if radius_min is not None and radius_min != Rmin:
raise ValueError(
"Rmin and radius_min have different values, "
"this is only allowed if one of them is None"
)
radius_min = Rmin
if radius_min is None:
raise ValueError("radius_min must be set")
if Rmax is not None:
warnings.warn(
"Rmax is deprecated, use radius_max instead",
category=FutureWarning,
)
if radius_max is not None and radius_max != Rmax:
raise ValueError(
"Rmax and radius_max have different values, "
"this is only allowed if one of them is None"
)
radius_max = Rmax
if radius_max is None:
raise ValueError("radius_max must be set")

self.targetN=targetN
self.convert_nbody=convert_nbody
self.densitypower=densitypower
self.thermalpower=thermalpower
self.Rmin=radius_min
self.Rmax=radius_max
self.gamma=gamma
self.q_out=q_out
self.discfraction=discfraction

self.a=self.thermalpower
self.a2=self.thermalpower/2
self.g=densitypower
self.g2=2-densitypower
self.k_out=((1+discfraction)/Rmax**3)**0.5
self.sigma_out=self.g2*discfraction/(2*numpy.pi*Rmax**self.g*(Rmax**self.g2-Rmin**self.g2))
self.cs_out=self.q_out*numpy.pi*self.sigma_out/self.k_out

self.base_cylinder=uniform_unit_cylinder(targetN,base_grid)


def sigma(self,r):
return self.sigma_out*(self.Rmax/r)**self.g

def csound(self,r):
return self.cs_out*(self.Rmax/r)**self.a2

def cmass(self,r):
return self.discfraction*(r**self.g2-self.Rmin**self.g2)/(self.Rmax**self.g2-self.Rmin**self.g2)

def mass_encl(self,r):
return 1+self.cmass(r)

def kappa(self,r):
return (self.mass_encl(r)/r**3)**0.5

def toomreQ(self,r):
return self.csound(r)*self.kappa(r)/numpy.pi/self.sigma(r)

def getradius(self,f):
return ((self.Rmax**self.g2-self.Rmin**self.g2)*f+self.Rmin**self.g2)**(1./self.g2)

def zscale(self,r):
return self.csound(r)/self.kappa(r)

def u(self,r):
if self.gamma ==1.:
self.targetN = targetN
self.convert_nbody = convert_nbody
self.densitypower = densitypower
self.thermalpower = thermalpower
self.Rmin = rename_fn_par("radius_min", radius_min, "Rmin", Rmin, 1)
self.Rmax = rename_fn_par("radius_max", radius_max, "Rmax", Rmax, 100)
self.gamma = gamma
self.q_out = q_out
self.discfraction = discfraction

self.a = self.thermalpower
self.a2 = self.thermalpower/2
self.g = densitypower
self.g2 = 2 - densitypower
self.k_out = ((1 + discfraction) / self.Rmax**3)**0.5
self.sigma_out = self.g2 * discfraction / (
2 * numpy.pi * self.Rmax**self.g *
(self.Rmax**self.g2 - self.Rmin**self.g2))
self.cs_out = self.q_out * numpy.pi * self.sigma_out / self.k_out

self.base_cylinder = uniform_unit_cylinder(targetN, base_grid)

@property
def radius_min(self):
return self.Rmin

@property
def radius_max(self):
return self.Rmax

def sigma(self, r):
return self.sigma_out * (self.Rmax / r)**self.g

def csound(self, r):
return self.cs_out * (self.Rmax / r)**self.a2

def cmass(self, r):
return self.discfraction * (r**self.g2 - self.Rmin**self.g2) / (
self.Rmax**self.g2 - self.Rmin**self.g2)

def mass_encl(self, r):
return 1 + self.cmass(r)

def kappa(self, r):
return (self.mass_encl(r) / r**3)**0.5

def toomreQ(self, r):
return self.csound(r) * self.kappa(r) / numpy.pi / self.sigma(r)

def getradius(self, f):
return (
(self.Rmax**self.g2 - self.Rmin**self.g2) * f + self.Rmin**self.g2)**(
1./self.g2)

def zscale(self, r):
return self.csound(r) / self.kappa(r)

def u(self, r):
if self.gamma == 1.:
return self.csound(r)**2
else:
return self.csound(r)**2/(self.gamma-1)
return self.csound(r)**2 / (self.gamma - 1)

def vcirc(self, r):
return (self.mass_encl(r) / r)**0.5

def vcirc(self,r):
return (self.mass_encl(r)/r)**0.5

def new_model(self):
x,y,z=self.base_cylinder.make_xyz()
self.actualN=len(x)
f=x**2+y**2
r=f**0.5
rtarget=self.getradius(f)
mass=self.discfraction*numpy.ones_like(x)/self.actualN
internal_energy=self.u(rtarget)
zscale=self.zscale(rtarget)
r=r.clip(1.e-8,2.)
x=x/r
y=y/r

vx=-y*self.vcirc(rtarget)
vy=x*self.vcirc(rtarget)
vz=numpy.zeros_like(x)
x=rtarget*x
y=rtarget*y
z=approximate_inverse_error_function(z)*zscale*2.**0.5
x, y, z = self.base_cylinder.make_xyz()
self.actualN = len(x)
f = x**2 + y**2
r = f**0.5
rtarget = self.getradius(f)

mass = self.discfraction * numpy.ones_like(x) / self.actualN
internal_energy = self.u(rtarget)
zscale = self.zscale(rtarget)
r = r.clip(1.e-8, 2.)
x= x / r
y= y / r

vx = -y * self.vcirc(rtarget)
vy = x * self.vcirc(rtarget)
vz = numpy.zeros_like(x)

x = rtarget * x
y = rtarget * y
z = approximate_inverse_error_function(z) * zscale * 2.**0.5

return (mass, x, y, z, vx, vy, vz, internal_energy)

@property
def result(self):
masses, x,y,z, vx,vy,vz, internal_energies = self.new_model()
masses, x, y, z, vx, vy, vz, internal_energies = self.new_model()
result = Particles(self.actualN)
result.mass = nbody_system.mass.new_quantity(masses)
result.x = nbody_system.length.new_quantity(x)
Expand All @@ -168,10 +160,10 @@ def result(self):
result.vy = nbody_system.speed.new_quantity(vy)
result.vz = nbody_system.speed.new_quantity(vz)
result.u = nbody_system.specific_energy.new_quantity(internal_energies)

if not self.convert_nbody is None:
result = ParticlesWithUnitsConverted(result, self.convert_nbody.as_converter_from_si_to_generic())
result = ParticlesWithUnitsConverted(
result, self.convert_nbody.as_converter_from_si_to_generic())
result = result.copy()

return result

return result
55 changes: 55 additions & 0 deletions src/amuse/support/helpers.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
import warnings


def rename_fn_par(new_name, new_value, old_name, old_value, default_value):
"""Get value of a renamed function parameter.

If you have a function f with parameter ``x`` and default value 1, like this:

.. code-block::

def f(x=1):
...

and you want to rename ``x`` to ``y`` while staying backwards compatible, then you
can rename ``x`` to ``y``, add ``x`` back in at the end so that old keyword
arguments still work, give both a default value of ``None``, and use this function
like this:

.. code-block::

def f(y=None, x=None):
value = rename_fn_par("y", y, "x", x, 1)

Any callers using ``x`` explicitly will receive a warning to change their code to
use ``y`` in the future. If both ``x`` and ``y`` are set and the values are
different, then an exception is raised.

Args:
new_name (str): The new name of the variable
new_value: The value passed using the new name
old_name (str): The old name of the variable
old_value: The value passed using the old name
default_value: The default value if neither are set

Returns:
Either new_value or old_value if only one is set or they're set to the same
value, or default_value if neither is set.

Raises:
ValueError: If both new_value and old_value are set, and to different values.
"""
if new_value is not None:
if old_value is not None and old_value != new_value:
raise ValueError(
f"{old_name} and {new_name} have different values,"
" which is not allowed because they represent the same thing.")
return new_value

if old_value is not None:
warnings.warn(
f"{old_name} is deprecated, please use {new_name} instead",
category=FutureWarning)
return old_value

return default_value
40 changes: 40 additions & 0 deletions src/tests/core_tests/test_helpers.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
from amuse.support.helpers import rename_fn_par

import pytest


@pytest.mark.filterwarnings("ignore: old is deprecated.*")
def test_rename_fn_par():
assert rename_fn_par("new", None, "old", None, 1) == 1
assert rename_fn_par("new", None, "old", 1, 2) == 1
assert rename_fn_par("new", 1, "old", None, 2) == 1
assert rename_fn_par("new", 1, "old", 1, 2) == 1

with pytest.raises(ValueError):
rename_fn_par("new", 1, "old", 2, 3)


class _RenamedMethod:
def __init__(self, new = None, a = 1, b = "test", old = None):
self.y = rename_fn_par("new", new, "old", old, 42)


@pytest.mark.filterwarnings("ignore: old is deprecated.*")
def test_rename_fn_par_usage():
rm = _RenamedMethod()
assert rm.y == 42

rm = _RenamedMethod(1)
assert rm.y == 1

rm = _RenamedMethod(new=1)
assert rm.y == 1

rm = _RenamedMethod(old=1)
assert rm.y == 1

rm = _RenamedMethod(1, 1, "test", 1)
assert rm.y == 1

with pytest.raises(ValueError):
_RenamedMethod(1, 2, "test", 3)
Loading