Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revive contact smoothing #989

Open
wants to merge 30 commits into
base: feature-contacts
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
b4786d5
hacky support for contact temperature smoothing
kecnry Oct 4, 2024
f0fe605
added .save_as_standard_mesh after smoothing
gecheline Nov 12, 2020
4720723
cache "smoothed_teffs"
kecnry Nov 13, 2020
a9c363f
update weighing to use highest teff2 instead of teff1_neck
gecheline Nov 13, 2020
7e429c7
added lateral, isotropic and spotty mixing
kecnry Oct 4, 2024
0bc0a0c
parameterization of mixing methods with user-provided power and fixed…
gecheline Nov 14, 2020
37a007d
updated limits for mixing power to None
gecheline Nov 14, 2020
1530d58
Added "spotted" option to mixing_method.
aprsa Nov 15, 2020
9b7116d
Commented out a print statement.
aprsa Nov 16, 2020
29c7a63
updates to distl and stuff to allow plotting
kecnry Oct 4, 2024
670df11
reinstate mixing code in backend
Nov 4, 2024
119a3ef
make function of two occurances of mixing in backends.py
Nov 6, 2024
3362d9a
rebuild default bundles
Nov 6, 2024
8ea7a0f
set default mix to false in component.py
Nov 6, 2024
7f2d846
implement perfect mixing
Nov 6, 2024
e8c4ff2
some documentation/description of the mixing methods
Nov 7, 2024
2b71ab1
update default bundles with new choice
Nov 7, 2024
aad340b
add simple test of perfect mixing
Nov 7, 2024
b97a6ea
Merge remote-tracking branch 'origin/feature-contacts' into revive-co…
matthiasfabry Dec 4, 2024
ae834bf
add documentation
matthiasfabry Dec 18, 2024
2b72d1b
pep8 formatting
matthiasfabry Dec 18, 2024
9a21b37
use logger; descriptive ValueErrors
matthiasfabry Dec 18, 2024
f321360
remove spotty transfer cause too parametric
matthiasfabry Dec 18, 2024
bfee696
gate energy transfer behind envelope requirement; rename to `energy_t…
matthiasfabry Dec 19, 2024
41c536c
add forbidden labels
matthiasfabry Dec 20, 2024
d227197
recompute def bundles
matthiasfabry Dec 20, 2024
75c2b5a
remove old implementation in universe.py
matthiasfabry Dec 20, 2024
d6a8d03
Revert "updates to distl and stuff to allow plotting"
matthiasfabry Dec 20, 2024
caa9fee
Revert "recompute def bundles"
matthiasfabry Feb 4, 2025
b425bb2
update only relevant parameters in default bundles
matthiasfabry Feb 4, 2025
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
Prev Previous commit
Next Next commit
remove old implementation in universe.py
  • Loading branch information
matthiasfabry committed Dec 20, 2024
commit 75c2b5abf8180f45390c77fc3ed547712445e14d
257 changes: 0 additions & 257 deletions phoebe/backend/universe.py
Original file line number Diff line number Diff line change
@@ -3331,260 +3331,3 @@ def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None):
return teffs

raise NotImplementedError("teffext=True not yet supported for pulsations")


# def dist(p1, p2):
# (x1, y1, z1), (x2, y2, z2) = p1, p2
# return ((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)**0.5

# class ContactSmoothing(Feature):
# def __init__(self, w=0.5, cutoff=0, **kwargs):
# self._w = w
# self._cutoff = cutoff

# @classmethod
# def from_bundle(cls, b, feature):
# """
# Initialize a Spot feature from the bundle.
# """

# # feature_ps = b.get_feature(feature=feature, **_skip_filter_checks)

# # colat = feature_ps.get_value(qualifier='colat', unit=u.rad, **_skip_filter_checks)
# # longitude = feature_ps.get_value(qualifier='long', unit=u.rad, **_skip_filter_checks)

# # if len(b.hierarchy.get_stars())>=2:
# # star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks)
# # orbit_ps = b.get_component(component=b.hierarchy.get_parent_of(feature_ps.component), **_skip_filter_checks)
# # # TODO: how should this handle dpdt?

# # # we won't use syncpar directly because that is defined wrt sidereal period and we want to make sure
# # # this translated to roche longitude correctly. In the non-apsidal motion case
# # # syncpar = period_anom_orb / period_star
# # period_anom_orb = orbit_ps.get_value(qualifier='period_anom', unit=u.d, **_skip_filter_checks)
# # period_star = star_ps.get_value(qualifier='period', unit=u.d, **_skip_filter_checks)
# # dlongdt = 2*pi * (period_anom_orb/period_star - 1) / period_anom_orb
# # else:
# # star_ps = b.get_component(component=feature_ps.component, **_skip_filter_checks)
# # dlongdt = star_ps.get_value(qualifier='freq', unit=u.rad/u.d, **_skip_filter_checks)
# # longitude += np.pi/2

# # radius = feature_ps.get_value(qualifier='radius', unit=u.rad, **_skip_filter_checks)
# # relteff = feature_ps.get_value(qualifier='relteff', unit=u.dimensionless_unscaled, **_skip_filter_checks)

# # t0 = b.get_value(qualifier='t0', context='system', unit=u.d, **_skip_filter_checks)

# # return cls(colat, longitude, dlongdt, radius, relteff, t0)

# @staticmethod
# def _isolate_neck(coords_all, teffs_all, cutoff = 0.,component=1, plot=False):
# if component == 1:
# cond = coords_all[:,0] >= 0+cutoff
# elif component == 2:
# cond = coords_all[:,0] <= 1-cutoff
# else:
# raise ValueError

# if plot:
# import matplotlib.pyplot as plt
# fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10,5))
# axes[0].scatter(coords_all[cond][:,0], coords_all[cond][:,1])
# axes[1].scatter(coords_all[cond][:,0], teffs_all[cond])
# axes[0].set_xlabel('x')
# axes[0].set_ylabel('y')
# axes[1].set_xlabel('x')
# axes[1].set_ylabel('teff')
# fig.tight_layout()
# plt.show()

# return coords_all[cond], teffs_all[cond]


# @staticmethod
# def _isolate_sigma_fitting_regions(coords_neck, teffs_neck, direction='x', cutoff=0, component=1, plot=False):
# from itertools import combinations

# distances = [dist(p1, p2) for p1, p2 in combinations(coords_neck, 2)]
# min_dist = np.min(distances[distances != 0])

# if direction == 'x':
# cond = (coords_neck[:,1] >= -0.2*min_dist) & (coords_neck[:,1] <= 0.2*min_dist)

# elif direction == 'y':

# if component == 1:
# cond = coords_neck[:,0] <= 0+cutoff+0.15*min_dist

# elif component == 2:
# cond = coords_neck[:,0] >= 1-cutoff-0.15*min_dist

# else:
# raise ValueError

# else:
# raise ValueError

# if plot:
# import matplotlib.pyplot as plt
# if direction == 'x':
# plt.scatter(coords_neck[cond][:,0], teffs_neck[cond])
# plt.show()
# elif direction == 'y':
# plt.scatter(coords_neck[cond][:,1], teffs_neck[cond])
# plt.show()


# return coords_neck[cond], teffs_neck[cond]

# @staticmethod
# def _compute_new_teff_at_neck(coords1, teffs1, coords2, teffs2, w=0.5, offset=0.):

# from itertools import combinations

# distances1 = [dist(p1, p2) for p1, p2 in combinations(coords1, 2)]
# min_dist1 = np.min(distances1[distances1 != 0])
# distances2 = [dist(p1, p2) for p1, p2 in combinations(coords2, 2)]
# min_dist2 = np.min(distances2[distances2 != 0])

# x_neck = np.average((coords1[:,0].max(), coords2[:,0].min()))
# amplitude_1 = teffs1.max()
# amplitude_2 = teffs2.max()

# teffs_neck1 = teffs1[coords1[:,0] >= coords1[:,0].max() - 0.25*min_dist1]
# teffs_neck2 = teffs2[coords2[:,0] <= coords2[:,0].min() + 0.25*min_dist2]

# teff1 = np.average(teffs_neck1)
# teff2 = np.average(teffs_neck2)
# tavg = w*teff1 + (1-w)*teff2
# if tavg > teffs2.max():
# print('Warning: Tavg > Teff2, setting new temperature to 1 percent of Teff2 max. %i > %i' % (int(tavg), int(teffs2.max())))
# tavg = teffs2.max() - 0.01*teffs2.max()

# return x_neck, tavg


# @staticmethod
# def _compute_sigmax(Tavg, x, x0, offset, amplitude):
# return ((-1)*(x-x0)**2/np.log((Tavg-offset)/amplitude))**0.5


# @staticmethod
# def _fit_sigma_amplitude(coords, teffs, offset, cutoff=0, direction='y', component=1, plot=False):

# def gaussian_1d(x, sigma):
# a = 1./sigma**2
# g = offset + amplitude * np.exp(- (a * ((x - x0) ** 2)))
# return g

# from scipy.optimize import curve_fit

# if direction == 'y':
# coord_ind = 1
# x0 = 0.
# elif direction == 'x':
# coord_ind = 0
# if component == 1:
# x0 = 0+cutoff
# elif component == 2:
# x0 = 1-cutoff
# else:
# raise ValueError
# else:
# raise ValueError

# amplitude = teffs.max() - offset
# sigma_0 = 0.5
# result = curve_fit(gaussian_1d,
# xdata=coords[:,coord_ind],
# ydata=teffs,
# p0=(0.5,),
# bounds=[0.01,1000])

# sigma = result[0]
# model = gaussian_1d(coords[:,coord_ind], sigma)

# if plot:
# import matplotlib.pyplot as plt
# plt.scatter(coords[:,coord_ind], teffs)
# plt.scatter(coords[:,coord_ind], model)
# plt.show()

# return sigma, amplitude, model


# @staticmethod
# def _compute_twoD_Gaussian(coords, sigma_x, sigma_y, amplitude, cutoff=0, offset=0, component=1):
# y0 = 0.
# x0 = 0+cutoff if component == 1 else 1-cutoff
# a = 1. / sigma_x ** 2
# b = 1. / sigma_y ** 2
# return offset + amplitude * np.exp(- (a * ((coords[:, 0] - x0) ** 2))) * np.exp(- (b * ((coords[:, 1] - y0) ** 2)))


# def smooth_teffs(self, xyz1, teffs1, xyz2, teffs2, offset=0.):

# coords_neck_1, teffs_neck_1 = self._isolate_neck(xyz1, teffs1, cutoff = self._cutoff, component=1, plot=False)
# coords_neck_2, teffs_neck_2 = self._isolate_neck(xyz2, teffs2, cutoff = self._cutoff, component=2, plot=False)

# x_neck, Tavg = self._compute_new_teff_at_neck(coords_neck_1, teffs_neck_1, coords_neck_2, teffs_neck_2, w=self._w, offset=offset)

# sigma_x1 = self._compute_sigmax(Tavg, x_neck, x0=0+self._cutoff, offset=offset, amplitude=teffs_neck_1.max())
# sigma_x2 = self._compute_sigmax(Tavg, x_neck, x0=1-self._cutoff, offset=offset, amplitude=teffs_neck_2.max())

# coords_fit_y1, teffs_fit_y1 = self._isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y', cutoff=self._cutoff,
# component=1, plot=False)
# coords_fit_y2, teffs_fit_y2 = self._isolate_sigma_fitting_regions(coords_neck_2, teffs_neck_2, direction='y', cutoff=self._cutoff,
# component=2, plot=False)

# sigma_y1, amplitude_y1, model_y1 = self._fit_sigma_amplitude(coords_fit_y1, teffs_fit_y1, offset, direction='y',
# component=1, plot=False)
# sigma_y2, amplitude_y2, model_y2 = self._fit_sigma_amplitude(coords_fit_y2, teffs_fit_y2, offset, direction='y',
# component=2, plot=False)

# new_teffs1 = self. _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(),
# cutoff=self._cutoff, offset=offset, component=1)
# new_teffs2 = self._compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(),
# cutoff=self._cutoff, offset=offset, component=2)

# return new_teffs1, new_teffs2


# @property
# def proto_coords(self):
# """
# """
# return True


# def process_teffs(self, teffs, coords, s=np.array([0., 0., 1.]), t=None):
# """
# Change the local effective temperatures for any values within the
# "cone" defined by the spot. Any teff within the spot will have its
# current value multiplied by the "relteff" factor

# :parameter array teffs: array of teffs for computations
# :parameter array coords: array of coords for computations
# :t float: current time
# """

# # THE COORDINATES NEED TO BE ROCHE!

# xyz1 = b['value@xyz_elements@mesh@primary'][:,0]
# teffs1 = b['value@teffs@mesh@primary']
# xyz2 = b['value@xyz_elements@mesh@secondary'][:,0]
# teffs2 = b['value@teffs@mesh@secondary']

# # if t is None:
# # # then assume at t0
# # t = self._t0

# # pointing_vector = self.pointing_vector(s,t)
# # logger.debug("spot.process_teffs at t={} with pointing_vector={} and radius={}".format(t, pointing_vector, self._radius))

# # cos_alpha_coords = np.dot(coords, pointing_vector) / np.linalg.norm(coords, axis=1)
# # cos_alpha_spot = np.cos(self._radius)

# # filter_ = cos_alpha_coords > cos_alpha_spot
# # teffs[filter_] = teffs[filter_] * self._relteff

# # return teffs