revive contact smoothing implementation

phoebe/backend/backends.py

@@ -14,7 +14,7 @@
 from phoebe.parameters import StringParameter, DictParameter, ArrayParameter, ParameterSet
 from phoebe.parameters.parameters import _extract_index_from_string
 from phoebe import dynamics
-from phoebe.backend import universe, horizon_analytic
+from phoebe.backend import universe, horizon_analytic, contacts_smoothing
 from phoebe.atmospheres import passbands
 from phoebe.distortions  import roche
 from phoebe.frontend import io
@@ -974,6 +974,31 @@ def _compute_intrinsic_system_at_t0(self, b, compute,
         # kinds = [b.get_dataset(dataset=ds).kind for ds in datasets]
 
         system.update_positions(t0, x0, y0, z0, vx0, vy0, vz0, etheta0, elongan0, eincl0, ignore_effects=True)
+
+        #TEMPERATURE SMOOTHING FOR CONTACTS
+        # smoothing_enabled = b.get_value(qualifier='smoothing_enabled', context='component', **_skip_filter_checks)
+
+        # if smoothing_enabled:
+        #     smoothing_method = b.get_value(qualifier='smoothing_method', context='component', **_skip_filter_checks)
+        #     smoothing_factor = b.get_value(qualifier='smoothing_factor', context='component', **_skip_filter_checks)
+        #     primary_mesh, secondary_mesh = system.bodies[0].meshes.values()
+        #     coords1 = primary_mesh.roche_coords_for_computations
+        #     teffs1 = primary_mesh.teffs
+        #     coords2 = secondary_mesh.roche_coords_for_computations
+        #     teffs2 = secondary_mesh.teffs
+
+        #     new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1),
+        #                                                             np.array(coords2), np.array(teffs2),
+        #                                                             smoothing_method=smoothing_method)
+        #                                                             # w=smoothing_factor, cutoff=0.)
+        #     # primary_mesh.update_columns(teffs=new_teffs1)
+        #     # secondary_mesh.update_columns(teffs=new_teffs2)
+
+        #     # print("new_teffs1.median", np.median(new_teffs1))
+        #     # print("new_teffs2.median", np.median(new_teffs2))
+        #     system.bodies[0]._halves[0].smoothed_teffs = new_teffs1
+        #     system.bodies[0]._halves[1].smoothed_teffs = new_teffs2
+
         system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True)
 
 
@@ -1070,6 +1095,36 @@ def _run_single_time(self, b, i, time, infolist, **kwargs):
             logger.debug("rank:{}/{} PhoebeBackend._run_single_time: calling system.update_positions at time={}".format(mpi.myrank, mpi.nprocs, time))
             system.update_positions(time, xi, yi, zi, vxi, vyi, vzi, ethetai, elongani, eincli, ds=di, Fs=Fi)
 
+            #TEMPERATURE SMOOTHING FOR CONTACTS
+            mixing_enabled = b.get_value(qualifier='mixing_enabled', context='component', **_skip_filter_checks)
+
+            if mixing_enabled and i==0:
+                mixing_method = b.get_value(qualifier='mixing_method', context='component', **_skip_filter_checks)
+                mixing_power = b.get_value(qualifier='mixing_power', context='component', **_skip_filter_checks)
+                secondary_teff = b.get_value(qualifier='teff', context='component', component='secondary', **_skip_filter_checks)
+                primary_teff = b.get_value(qualifier='teff', context='component', component='primary', **_skip_filter_checks)
+                teff_factor = 1.-secondary_teff/primary_teff
+
+                primary_mesh, secondary_mesh = system.bodies[0].meshes.values()
+                coords1 = primary_mesh.roche_coords_for_computations
+                teffs1 = primary_mesh.teffs
+                coords2 = secondary_mesh.roche_coords_for_computations
+                teffs2 = secondary_mesh.teffs
+
+                new_teffs1, new_teffs2 = contacts_smoothing.smooth_teffs(np.array(coords1), np.array(teffs1),
+                                                                        np.array(coords2), np.array(teffs2),
+                                                                        mixing_method=mixing_method,
+                                                                        mixing_power=mixing_power, teff_factor=teff_factor)
+                                                                        # w=smoothing_factor, cutoff=0.)
+                # primary_mesh.update_columns(teffs=new_teffs1)
+                # secondary_mesh.update_columns(teffs=new_teffs2)
+
+                # print("new_teffs1.median", np.median(new_teffs1))
+                # print("new_teffs2.median", np.median(new_teffs2))
+                system.bodies[0]._halves[0].smoothed_teffs = new_teffs1
+                system.bodies[0]._halves[1].smoothed_teffs = new_teffs2
+
+
             # Now we need to determine which triangles are visible and handle subdivision
             # NOTE: this should come after populate_observables so that each subdivided triangle
             # will have identical local quantities.  The only downside to this is that we can't

phoebe/backend/contacts_smoothing.py

@@ -0,0 +1,227 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from itertools import combinations
+
+
+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:
+        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], np.argwhere(cond).flatten()
+
+def _dist(p1, p2):
+    (x1, y1, z1), (x2, y2, z2) = p1, p2
+    return ((x2 - x1)**2 + (y2 - y1)**2 + (z2 - z1)**2)**0.5
+
+
+def _isolate_sigma_fitting_regions(coords_neck, teffs_neck, direction='x', cutoff=0., component=1, plot=False):
+
+    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:
+        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]
+
+def _compute_new_teff_at_neck(coords1, teffs1, coords2, teffs2, w=0.5, offset=0.):
+
+    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.max(teffs2)#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
+
+
+def _compute_sigmax(Tavg, x, x0, offset, amplitude):
+    return ((-1)*(x-x0)**2/np.log((Tavg-offset)/amplitude))**0.5
+
+
+def _fit_sigma_amplitude(coords, teffs, offset=0., 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:
+        plt.scatter(coords[:,coord_ind], teffs)
+        plt.scatter(coords[:,coord_ind], model)
+        plt.show()
+
+    return sigma, amplitude, model
+
+
+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 gaussian_smoothing(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.):
+    coords_neck_1, teffs_neck_1, cond1 = _isolate_neck(xyz1, teffs1, cutoff = cutoff, component=1, plot=False)
+    coords_neck_2, teffs_neck_2, cond2 = _isolate_neck(xyz2, teffs2, cutoff = cutoff, component=2, plot=False)
+
+    x_neck, Tavg = _compute_new_teff_at_neck(coords_neck_1, teffs_neck_1, coords_neck_2, teffs_neck_2, w=w, offset=offset)
+
+    sigma_x1 = _compute_sigmax(Tavg, x_neck, x0=0+cutoff, offset=offset, amplitude=teffs_neck_1.max())
+    sigma_x2 = _compute_sigmax(Tavg, x_neck, x0=1-cutoff, offset=offset, amplitude=teffs_neck_2.max())
+
+    coords_fit_y1, teffs_fit_y1 = _isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y', cutoff=cutoff, 
+                                                                        component=1, plot=False)
+    coords_fit_y2, teffs_fit_y2 = _isolate_sigma_fitting_regions(coords_neck_2, teffs_neck_2, direction='y', cutoff=cutoff,
+                                                                        component=2, plot=False)
+
+    sigma_y1, amplitude_y1, model_y1 = _fit_sigma_amplitude(coords_fit_y1, teffs_fit_y1, offset, direction='y', 
+                                                                        component=1, plot=False)
+    sigma_y2, amplitude_y2, model_y2 = _fit_sigma_amplitude(coords_fit_y2, teffs_fit_y2, offset, direction='y', 
+                                                                        component=2, plot=False)
+
+    new_teffs1 =  _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(), 
+                                                cutoff=cutoff, offset=offset, component=1)
+    new_teffs2 = _compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(), 
+                                                cutoff=cutoff, offset=offset, component=2)
+
+    # print(cond1, len(cond1), len(new_teffs1))
+    # print(cond2, len(cond2), len(new_teffs2))
+    teffs1[cond1] = new_teffs1
+    teffs2[cond2] = new_teffs2
+
+    return teffs1, teffs2
+
+
+def lateral_transfer(t2s, teffs2, mixing_power=0.5, teff_factor=1.):
+    x2s = t2s[:,0]
+    y2s = t2s[:,1]
+    z2s = t2s[:,2]
+
+    y2s_neck = y2s[x2s<1]
+    z2s_neck = z2s[x2s<1]
+    rs_neck = (y2s_neck**2+z2s_neck**2)**0.5
+    lat = np.min(rs_neck)
+    # print('mixing latitude = ', lat)
+    # lat = 0.15
+    # power = 0.75
+    filt = (z2s>-lat) & (z2s<lat)
+    c = (lat-np.abs(z2s[filt]))**mixing_power
+    teffs2[filt] *= 1+teff_factor*c/c.max()  
+
+    return teffs2
+
+def isotropic_transfer(t2s, teffs2, mixing_power=0.5, teff_factor=1.):
+    d2s = np.sqrt(t2s[:,0]*t2s[:,0] + t2s[:,1]*t2s[:,1] + t2s[:,2]*t2s[:,2])
+    teffs2 *= 1+teff_factor*(1-((d2s-d2s.min())/(d2s.max()-d2s.min()))**mixing_power)
+
+    return teffs2
+
+
+def spotty_transfer(t2s, teffs2):
+    d2s = np.sqrt(t2s[:,0]*t2s[:,0] + t2s[:,1]*t2s[:,1] + t2s[:,2]*t2s[:,2])
+    for s in range(10):
+        idx = int(len(d2s)*np.random.rand())
+        size = 0.3*np.random.rand()
+        factor = 1.1-0.025*np.random.rand()
+
+        ds = np.sqrt((t2s[:,0]-t2s[idx,0])**2 + (t2s[:,1]-t2s[idx,1])**2 + (t2s[:,2]-t2s[idx,2])**2)
+        teffs2[ds<size] *= factor
+
+    return teffs2
+
+
+def smooth_teffs(xyz1, teffs1, xyz2, teffs2, mixing_method='lateral', mixing_power=0.5, teff_factor=1.):
+
+    if mixing_method == 'lateral':
+        teffs2 = lateral_transfer(xyz2, teffs2, mixing_power, teff_factor)
+
+    elif mixing_method == 'isotropic':
+        teffs2 = isotropic_transfer(xyz2, teffs2, mixing_power, teff_factor)
+
+    elif mixing_method == 'spotty':
+        teffs2 = spotty_transfer(xyz2, teffs2)
+
+    else:
+        teffs1, teffs2 = gaussian_smoothing(xyz1, teffs1, xyz2, teffs2)
+
+    return teffs1, teffs2