Skip to content

Commit ff24a8d

Browse files
spzwartspzwart
committed
Added high order bridge
1 parent 96ba87b commit ff24a8d

1 file changed

Lines changed: 243 additions & 0 deletions

File tree

‎high_order_bridge.ipynb‎

Lines changed: 243 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,243 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"AMUSE tutorial on high-order bridge\n",
8+
"====================\n",
9+
"\n",
10+
"Hierarchical coupling strategies are fundamental parts of AMUSE.\n",
11+
"It enables us to combined the output of a wide variety of solvers into a homogeneous solution.\n",
12+
"In this example we will be nesting multiple bridges, to show the power of bridge."
13+
]
14+
},
15+
{
16+
"cell_type": "code",
17+
"execution_count": null,
18+
"metadata": {},
19+
"outputs": [],
20+
"source": [
21+
"import numpy\n",
22+
"from amuse.units import (units, constants)\n",
23+
"#from amuse.lab import Particles"
24+
]
25+
},
26+
{
27+
"cell_type": "code",
28+
"execution_count": null,
29+
"metadata": {},
30+
"outputs": [],
31+
"source": [
32+
"def orbital_period(Mtot, a):\n",
33+
" return (((4 * numpy.pi**2) * a**3)/(constants.G * Mtot)).sqrt()\n",
34+
"\n",
35+
"from amuse.ext.orbital_elements import new_binary_from_orbital_elements\n",
36+
"def get_star_planet_and_moon():\n",
37+
" Mstar = 1.0|units.MSun\n",
38+
" Mplanet = 1.0|units.MJupiter\n",
39+
" a_planet = 5.2 | units.au\n",
40+
" e_planet = 0.0\n",
41+
" bodies = new_binary_from_orbital_elements(Mstar, Mplanet, a_planet, e_planet,\n",
42+
" G=constants.G)\n",
43+
" bodies[0].name = \"star\"\n",
44+
" bodies[1].name = \"planet\"\n",
45+
" planet = bodies[bodies.name==\"planet\"]\n",
46+
" RH_planet = a_planet*(1.0-e_planet)*(Mplanet/(3*Mstar))**(1./3.)\n",
47+
" a_moon = 0.1*RH_planet\n",
48+
" e_moon = 0.0\n",
49+
" Mmoon = 0.01*Mplanet\n",
50+
" Pmoon = orbital_period(Mplanet+Mmoon, a_moon)\n",
51+
" moon = new_binary_from_orbital_elements(planet.mass, \n",
52+
" Mmoon, a_moon, e_moon,\n",
53+
" G=constants.G)\n",
54+
" moon.position -= moon[0].position\n",
55+
" moon.velocity -= moon[0].velocity\n",
56+
" moon = moon[1].as_set()\n",
57+
" moon.position += planet.position\n",
58+
" moon.velocity += planet.velocity\n",
59+
" moon.semimajor_axis = a_moon\n",
60+
" moon.name = \"moon\"\n",
61+
" bodies.add_particle(moon)\n",
62+
" bodies.move_to_center()\n",
63+
" return bodies\n",
64+
"bodies = get_star_planet_and_moon()\n",
65+
"print(bodies)"
66+
]
67+
},
68+
{
69+
"cell_type": "markdown",
70+
"metadata": {},
71+
"source": [
72+
"Now we have the orbits of the three particles, star, planet and moon. We now want to make a disk around the moon."
73+
]
74+
},
75+
{
76+
"cell_type": "code",
77+
"execution_count": null,
78+
"metadata": {},
79+
"outputs": [],
80+
"source": [
81+
"from amuse.ext.protodisk import ProtoPlanetaryDisk\n",
82+
"\n",
83+
"from amuse.lab import nbody_system\n",
84+
"\n",
85+
"converter = nbody_system.nbody_to_si(moon.mass.sum(), \n",
86+
" 1|units.au)\n",
87+
"def make_disk_around_celestial_body(moon, Mplanet):\n",
88+
"\n",
89+
" R = 1|units.au\n",
90+
" a_moon = moon.semimajor_axis\n",
91+
" e_moon = 0.0\n",
92+
" Mmoon = moon.mass.sum()\n",
93+
" RH_moon = a_moon*(1.0-e_moon)*(Mmoon/(3*Mplanet))**(1./3.)\n",
94+
" converter = nbody_system.nbody_to_si(Mmoon, R)\n",
95+
" Ndisk = 1000\n",
96+
" Rin = 0.03*RH_moon\n",
97+
" Rout = 0.3*RH_moon\n",
98+
" Pinner = orbital_period(Mmoon, Rin)\n",
99+
" Mdisk = 0.01 * Mmoon\n",
100+
"\n",
101+
" disk = ProtoPlanetaryDisk(Ndisk,\n",
102+
" convert_nbody=converter,\n",
103+
" Rmin=Rin/R,\n",
104+
" Rmax=Rout/R,\n",
105+
" q_out=10.0,\n",
106+
" discfraction=Mdisk/Mmoon).result\n",
107+
" disk.move_to_center()\n",
108+
" disk.position += moon.position\n",
109+
" disk.velocity += moon.velocity\n",
110+
"\n",
111+
" masses = Mdisk/float(Ndisk)\n",
112+
" disk.mass = masses\n",
113+
" rho = 3.0| (units.g/units.cm**3)\n",
114+
" disk.radius = (disk.mass/(4*rho))**(1./3.)\n",
115+
" return disk, Pinner\n",
116+
"\n",
117+
"planet = bodies[bodies.name==\"planet\"]\n",
118+
"moon = bodies[bodies.name==\"moon\"]\n",
119+
"disk, Pinner = make_disk_around_celestial_body(moon, \n",
120+
" planet.mass.sum())"
121+
]
122+
},
123+
{
124+
"cell_type": "code",
125+
"execution_count": null,
126+
"metadata": {},
127+
"outputs": [],
128+
"source": [
129+
"from amuse.community.huayno.interface import Huayno\n",
130+
"gravityA = Huayno(converter)\n",
131+
"gravityA.particles.add_particles(bodies)\n",
132+
"channel = {\"from stars\": bodies.new_channel_to(gravityA.particles),\n",
133+
" \"to_stars\": gravityA.particles.new_channel_to(bodies)}\n",
134+
"\n",
135+
"gravityB = Huayno(converter, mode=\"openmp\")\n",
136+
"gravityB.particles.add_particles(disk)\n",
137+
"channel.update({\"from_disk\": disk.new_channel_to(gravityB.particles)})\n",
138+
"channel.update({\"to_disk:\": gravityB.particles.new_channel_to(disk)})\n",
139+
"bodies.add_particles(disk)"
140+
]
141+
},
142+
{
143+
"cell_type": "code",
144+
"execution_count": null,
145+
"metadata": {},
146+
"outputs": [],
147+
"source": [
148+
"from amuse.couple import bridge\n",
149+
"from amuse.ext.composition_methods import *\n",
150+
"gravhydro = bridge.Bridge(use_threading=False, method=SPLIT_4TH_S_M4)\n",
151+
"gravhydro.add_system(gravityA, (gravityB,))\n",
152+
"gravhydro.add_system(gravityB, (gravityA,))\n",
153+
"gravhydro.timestep = 0.5*Pinner"
154+
]
155+
},
156+
{
157+
"cell_type": "code",
158+
"execution_count": null,
159+
"metadata": {},
160+
"outputs": [],
161+
"source": [
162+
"from amuse.ext.composition_methods import *\n",
163+
"from amuse.ext.orbital_elements import orbital_elements_from_binary\n",
164+
"\n",
165+
"def gravity_hydro_bridge(gravityA, gravityB, gravhydro, bodies,\n",
166+
" t_end):\n",
167+
"\n",
168+
" gravity_initial_total_energy = gravityA.get_total_energy() + gravityB.get_total_energy()\n",
169+
" model_time = 0 | units.Myr\n",
170+
" dt = 0.012|units.yr #1.0*Pinner\n",
171+
" while model_time < t_end:\n",
172+
"\n",
173+
" model_time += dt\n",
174+
"\n",
175+
" orbit_planet = orbital_elements_from_binary(bodies[:2], G=constants.G)\n",
176+
" orbit_moon = orbital_elements_from_binary(bodies[1:3], G=constants.G)\n",
177+
" print(\"Planet:\", \"ae=\", orbit_planet[2].in_(units.AU), orbit_planet[3])\n",
178+
" print(\"Moon:\", \"ae=\", orbit_moon[2].in_(units.AU), orbit_moon[3])\n",
179+
" \n",
180+
" dE_gravity = gravity_initial_total_energy/(gravityA.get_total_energy()+gravityB.get_total_energy())\n",
181+
" print(\"Time:\", model_time.in_(units.day), \\\n",
182+
" \"dE=\", dE_gravity)#, dE_hydro\n",
183+
"\n",
184+
" gravhydro.evolve_model(model_time)\n",
185+
" channel[\"to_stars\"].copy()\n",
186+
" channel[\"to_disk\"].copy()\n",
187+
" print(\"S=\", bodies[:3])\n",
188+
" print(\"g=\", gravityA.particles)\n",
189+
" print(gravityA.particles.y.in_(units.au), stars.y.in_(units.au))\n",
190+
" \n",
191+
" gravityA.stop()\n",
192+
" gravityB.stop()\n",
193+
"\n",
194+
"t_end = 1.0 | units.yr\n",
195+
"gravity_hydro_bridge(gravityA, gravityB, gravhydro, \n",
196+
" bodies, t_end)"
197+
]
198+
},
199+
{
200+
"cell_type": "markdown",
201+
"metadata": {},
202+
"source": [
203+
"You have created a \n",
204+
"\n",
205+
"Assignmnets and questions:\n",
206+
"---------------\n",
207+
"\n",
208+
"### Assignment 1:\n",
209+
"\n",
210+
"\n",
211+
"### Question 1:\n"
212+
]
213+
},
214+
{
215+
"cell_type": "code",
216+
"execution_count": null,
217+
"metadata": {},
218+
"outputs": [],
219+
"source": []
220+
}
221+
],
222+
"metadata": {
223+
"kernelspec": {
224+
"display_name": "Python 3",
225+
"language": "python",
226+
"name": "python3"
227+
},
228+
"language_info": {
229+
"codemirror_mode": {
230+
"name": "ipython",
231+
"version": 3
232+
},
233+
"file_extension": ".py",
234+
"mimetype": "text/x-python",
235+
"name": "python",
236+
"nbconvert_exporter": "python",
237+
"pygments_lexer": "ipython3",
238+
"version": "3.7.6"
239+
}
240+
},
241+
"nbformat": 4,
242+
"nbformat_minor": 4
243+
}

0 commit comments

Comments
 (0)