0
# Grid and Stellar Functions
1
2
Functions for setting up computational grids and stellar properties in disk simulations. These functions provide essential infrastructure for the spatial discretization and stellar parameter calculations in DustPy simulations.
3
4
## Capabilities
5
6
### Grid Functions
7
8
Functions for managing computational grids and derived grid quantities.
9
10
```python { .api }
11
from dustpy.std import grid
12
13
def grid.OmegaK(sim):
14
"""
15
Calculate Keplerian orbital frequency.
16
17
Computes the Keplerian frequency at each radial grid point
18
based on the stellar mass and radial position. This is
19
fundamental for calculating orbital time scales and
20
gas pressure support.
21
22
Parameters:
23
- sim: Simulation object
24
25
Returns:
26
numpy.ndarray: Keplerian frequency [1/s] at each radial point
27
28
Formula:
29
Ω_K = √(GM_star / r³)
30
31
Where:
32
- G: Gravitational constant
33
- M_star: Stellar mass
34
- r: Radial distance from star
35
"""
36
```
37
38
### Stellar Functions
39
40
Functions for calculating stellar properties and luminosity.
41
42
```python { .api }
43
from dustpy.std import star
44
45
def star.luminosity(sim):
46
"""
47
Calculate stellar luminosity.
48
49
Computes the stellar luminosity based on the stellar
50
mass, radius, and effective temperature using the
51
Stefan-Boltzmann law.
52
53
Parameters:
54
- sim: Simulation object
55
56
Returns:
57
float: Stellar luminosity [erg/s]
58
59
Formula:
60
L = 4π R²_star σ_SB T⁴_eff
61
62
Where:
63
- R_star: Stellar radius
64
- σ_SB: Stefan-Boltzmann constant
65
- T_eff: Effective temperature
66
"""
67
```
68
69
## Usage Examples
70
71
### Basic Grid Setup
72
73
```python
74
from dustpy import Simulation
75
from dustpy.std import grid
76
import dustpy.constants as c
77
78
# Create simulation with custom grid
79
sim = Simulation()
80
81
# Set grid parameters
82
sim.ini.grid.Nr = 200 # Number of radial points
83
sim.ini.grid.rmin = 0.1 * c.au # Inner radius
84
sim.ini.grid.rmax = 1000 * c.au # Outer radius
85
86
# Create grids
87
sim.makegrids()
88
89
# Calculate Keplerian frequency
90
omega_k = grid.OmegaK(sim) # [1/s]
91
92
# Convert to useful units
93
orbital_periods = 2 * c.pi / omega_k # [s]
94
periods_years = orbital_periods / c.year # [years]
95
96
print(f"Orbital period at 1 AU: {periods_years[50]:.1f} years")
97
print(f"Orbital period at 10 AU: {periods_years[150]:.1f} years")
98
```
99
100
### Stellar Property Calculations
101
102
```python
103
from dustpy.std import star
104
import dustpy.constants as c
105
106
# Set stellar parameters
107
sim.ini.star.M = 0.5 * c.M_sun # Half solar mass
108
sim.ini.star.R = 0.7 * c.R_sun # 70% solar radius
109
sim.ini.star.T = 4500 # Effective temperature [K]
110
111
# Initialize simulation
112
sim.makegrids()
113
sim.initialize()
114
115
# Calculate luminosity
116
luminosity = star.luminosity(sim) # [erg/s]
117
solar_luminosity = 3.828e33 # [erg/s]
118
L_ratio = luminosity / solar_luminosity
119
120
print(f"Stellar luminosity: {luminosity:.2e} erg/s")
121
print(f"Luminosity ratio (L/L_sun): {L_ratio:.2f}")
122
```
123
124
### Grid-Dependent Calculations
125
126
```python
127
# Access grid properties
128
radial_centers = sim.grid.r # [cm]
129
radial_interfaces = sim.grid.ri # [cm]
130
annulus_areas = sim.grid.A # [cm²]
131
132
# Calculate orbital frequencies
133
omega_k = grid.OmegaK(sim) # [1/s]
134
135
# Derived time scales
136
dynamical_time = 1.0 / omega_k # [s]
137
hill_sphere_radius = radial_centers * (sim.star.M / (3 * sim.star.M))**(1/3)
138
139
# Print grid information
140
print(f"Grid extends from {sim.grid.r[0]/c.au:.2f} to {sim.grid.r[-1]/c.au:.1f} AU")
141
print(f"Number of grid points: {len(sim.grid.r)}")
142
print(f"Innermost dynamical time: {dynamical_time[0]/c.year:.2e} years")
143
```
144
145
### Stellar Evolution Effects
146
147
```python
148
# Time-dependent stellar properties
149
def evolve_stellar_properties(sim, time):
150
"""Update stellar properties as function of time."""
151
152
# Example: Simple stellar evolution model
153
age_myr = time / (1e6 * c.year) # Age in Myr
154
155
# Luminosity decreases with time for low-mass stars
156
if age_myr < 10:
157
L_factor = 1.5 * (1 + np.exp(-age_myr/2))
158
else:
159
L_factor = 1.0
160
161
# Update stellar luminosity
162
base_luminosity = 4 * c.pi * sim.star.R**2 * c.sigma_sb * sim.star.T**4
163
sim.star.L = base_luminosity * L_factor
164
165
return sim.star.L
166
167
# Apply stellar evolution
168
current_time = 5e6 * c.year # 5 Myr
169
new_luminosity = evolve_stellar_properties(sim, current_time)
170
print(f"Evolved luminosity: {new_luminosity:.2e} erg/s")
171
```
172
173
### Grid Resolution Analysis
174
175
```python
176
def analyze_grid_resolution(sim):
177
"""Analyze grid resolution and recommend improvements."""
178
179
r = sim.grid.r # [cm]
180
omega_k = grid.OmegaK(sim) # [1/s]
181
182
# Grid spacing
183
dr = np.diff(sim.grid.ri) # [cm]
184
relative_spacing = dr[:-1] / r[1:] # Relative grid spacing
185
186
# Hill sphere resolution
187
hill_radius = r * (sim.star.M / (3 * sim.star.M))**(1/3)
188
points_per_hill = hill_radius / dr
189
190
# Report
191
print(f"Average relative spacing: {np.mean(relative_spacing):.3f}")
192
print(f"Maximum relative spacing: {np.max(relative_spacing):.3f}")
193
print(f"Minimum Hill sphere resolution: {np.min(points_per_hill):.1f} points")
194
195
# Recommendations
196
if np.max(relative_spacing) > 0.1:
197
print("Warning: Grid spacing may be too coarse")
198
if np.min(points_per_hill) < 5:
199
print("Warning: Hill sphere under-resolved")
200
201
# Check grid quality
202
analyze_grid_resolution(sim)
203
```
204
205
### Custom Grid Functions
206
207
```python
208
def calculate_epicyclic_frequency(sim):
209
"""Calculate epicyclic frequency for eccentric orbits."""
210
211
# For Keplerian potential, epicyclic frequency equals Keplerian
212
omega_k = grid.OmegaK(sim)
213
kappa = omega_k # Epicyclic frequency
214
215
return kappa
216
217
def calculate_vertical_frequency(sim):
218
"""Calculate vertical oscillation frequency."""
219
220
# For thin disk approximation
221
omega_k = grid.OmegaK(sim)
222
nu_z = omega_k # Vertical frequency ≈ Keplerian frequency
223
224
return nu_z
225
226
# Calculate additional frequencies
227
epicyclic_freq = calculate_epicyclic_frequency(sim)
228
vertical_freq = calculate_vertical_frequency(sim)
229
230
# Oscillation periods
231
epicyclic_period = 2 * c.pi / epicyclic_freq
232
vertical_period = 2 * c.pi / vertical_freq
233
```
234
235
### Integration with Disk Physics
236
237
```python
238
# Combine grid and stellar functions with disk physics
239
from dustpy.std import gas
240
241
# Initialize simulation
242
sim.makegrids()
243
sim.initialize()
244
245
# Calculate grid-dependent quantities
246
omega_k = grid.OmegaK(sim) # [1/s]
247
luminosity = star.luminosity(sim) # [erg/s]
248
249
# Use in disk temperature calculation
250
temperature = (luminosity / (16 * c.pi * c.sigma_sb * sim.grid.r**2))**0.25
251
sim.gas.T[:] = temperature
252
253
# Calculate pressure scale height
254
cs = gas.cs_isothermal(sim) # [cm/s]
255
scale_height = cs / omega_k # [cm]
256
257
print(f"Temperature at 1 AU: {temperature[50]:.0f} K")
258
print(f"Scale height at 1 AU: {scale_height[50]/c.au:.3f} AU")
259
```
260
261
### Grid Diagnostics
262
263
```python
264
def grid_diagnostics(sim):
265
"""Comprehensive grid quality diagnostics."""
266
267
r = sim.grid.r
268
ri = sim.grid.ri
269
270
print("Grid Diagnostics:")
271
print(f" Radial range: {r[0]/c.au:.2f} - {r[-1]/c.au:.0f} AU")
272
print(f" Number of points: {len(r)}")
273
274
# Grid spacing
275
dr = np.diff(ri)
276
print(f" Grid spacing: {dr[0]/c.au:.4f} - {dr[-1]/c.au:.2f} AU")
277
278
# Logarithmic spacing check
279
r_ratio = r[1:] / r[:-1]
280
is_logarithmic = np.allclose(r_ratio, r_ratio[0], rtol=0.01)
281
print(f" Logarithmic spacing: {is_logarithmic}")
282
283
if is_logarithmic:
284
print(f" Log spacing factor: {r_ratio[0]:.3f}")
285
286
# Time scale resolution
287
omega_k = grid.OmegaK(sim)
288
min_period = 2 * c.pi / np.max(omega_k) / c.year
289
max_period = 2 * c.pi / np.min(omega_k) / c.year
290
291
print(f" Orbital periods: {min_period:.1f} - {max_period:.0f} years")
292
293
# Run diagnostics
294
grid_diagnostics(sim)
295
```