From 9480b851b005671c92970774e718d6f515f7b8b4 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 11 Jan 2019 14:50:03 +0100
Subject: [PATCH] New phase field step that uses 4th order FD to compute force
- complicated pressure tensor derivation not required
---
phasefield/n_phase_boyer.py | 18 +--
phasefield/phasefieldstep_direct.py | 175 ++++++++++++++++++++++++++++
2 files changed, 181 insertions(+), 12 deletions(-)
create mode 100644 phasefield/phasefieldstep_direct.py
diff --git a/phasefield/n_phase_boyer.py b/phasefield/n_phase_boyer.py
index a1e2eec7..75f88f50 100644
--- a/phasefield/n_phase_boyer.py
+++ b/phasefield/n_phase_boyer.py
@@ -173,18 +173,12 @@ class capital_h(sp.Function):
else:
raise sp.function.ArgumentIndexError(self, argindex)
- @staticmethod
- def insert(expr):
- if isinstance(expr, capital_h):
- u, v = expr.args
- av = sp.Abs(v)
- zero_cond = sp.And(sp.Eq(u, 0), sp.Eq(v, 0))
- return sp.Piecewise((0, zero_cond),
- (av * v / (av + u**2), True))
- else:
- new_args = [capital_h.insert(arg) for arg in expr.args]
- return expr.func(*new_args) if new_args else expr
-
+ def doit(self, **hints):
+ u, v = self.args
+ av = sp.Abs(v)
+ zero_cond = sp.And(sp.Eq(u, 0), sp.Eq(v, 0))
+ return sp.Piecewise((0, zero_cond),
+ (av * v / (av + u**2), True))
def correction_g(c, surface_tensions, symbolic_coefficients=False):
assert len(c) == surface_tensions.rows
diff --git a/phasefield/phasefieldstep_direct.py b/phasefield/phasefieldstep_direct.py
new file mode 100644
index 00000000..29f3a196
--- /dev/null
+++ b/phasefield/phasefieldstep_direct.py
@@ -0,0 +1,175 @@
+import numpy as np
+from lbmpy.creationfunctions import create_lb_function
+from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
+from lbmpy.phasefield.analytical import force_from_phi_and_mu
+from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
+from lbmpy.phasefield.kerneleqs import mu_kernel
+from lbmpy.stencils import get_stencil
+from pystencils import create_data_handling, Assignment, create_kernel
+from pystencils.fd import discretize_spatial
+from pystencils.fd.spatial import fd_stencils_forth_order_isotropic
+from pystencils.simp import sympy_cse_on_assignment_list
+from pystencils.slicing import SlicedGetterDataHandling
+
+
+class PhaseFieldStepDirect:
+
+ def __init__(self, free_energy, order_parameters, domain_size, data_handling=None, name='pfn',
+ hydro_dynamic_relaxation_rate=1.0,
+ cahn_hilliard_relaxation_rates=1.0,
+ cahn_hilliard_gammas=1,
+ optimization=None):
+ if optimization is None:
+ optimization = {'openmp': False, 'target': 'cpu'}
+ openmp = optimization.get('openmp', False)
+ target = optimization.get('target', 'cpu')
+
+ if not hasattr(cahn_hilliard_relaxation_rates, '__len__'):
+ cahn_hilliard_relaxation_rates = [cahn_hilliard_relaxation_rates] * len(order_parameters)
+
+ if not hasattr(cahn_hilliard_gammas, '__len__'):
+ cahn_hilliard_gammas = [cahn_hilliard_gammas] * len(order_parameters)
+
+ if data_handling is None:
+ data_handling = create_data_handling(domain_size, periodicity=True, parallel=False)
+ dh = data_handling
+ self.data_handling = dh
+
+ stencil = get_stencil('D3Q19' if dh.dim == 3 else 'D2Q9')
+
+ self.free_energy = free_energy
+
+ # Data Handling
+ kernel_parameters = {'cpu_openmp': openmp, 'target': target, 'ghost_layers': 2}
+ gpu = target == 'gpu'
+ phi_size = len(order_parameters)
+ gl = kernel_parameters['ghost_layers']
+ self.phi_field = dh.add_array('{}_phi'.format(name), values_per_cell=phi_size, gpu=gpu, latex_name='φ',
+ ghost_layers=gl)
+ self.mu_field = dh.add_array("{}_mu".format(name), values_per_cell=phi_size, gpu=gpu, latex_name="μ",
+ ghost_layers=gl)
+ self.vel_field = dh.add_array("{}_u".format(name), values_per_cell=data_handling.dim, gpu=gpu, latex_name="u",
+ ghost_layers=gl)
+
+ self.force_field = dh.add_array("{}_force".format(name), values_per_cell=dh.dim, gpu=gpu, latex_name="F",
+ ghost_layers=gl)
+
+ self.phi = SlicedGetterDataHandling(self.data_handling, self.phi_field.name)
+ self.mu = SlicedGetterDataHandling(self.data_handling, self.mu_field.name)
+ self.velocity = SlicedGetterDataHandling(self.data_handling, self.vel_field.name)
+ self.force = SlicedGetterDataHandling(self.data_handling, self.force_field.name)
+
+ self.ch_pdfs = []
+ for i in range(len(order_parameters)):
+ src = dh.add_array("{}_ch_src{}".format(name, i), values_per_cell=len(stencil), ghost_layers=gl)
+ dst = dh.add_array("{}_ch_dst{}".format(name, i), values_per_cell=len(stencil), ghost_layers=gl)
+ self.ch_pdfs.append((src, dst))
+ self.hydro_pdfs = (dh.add_array("{}_hydro_src".format(name), values_per_cell=len(stencil), ghost_layers=gl),
+ dh.add_array("{}_hydro_dst".format(name), values_per_cell=len(stencil), ghost_layers=gl))
+
+ # Compute Kernels
+ mu_assignments = mu_kernel(self.free_energy, order_parameters, self.phi_field, self.mu_field)
+ mu_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in mu_assignments]
+ mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
+ self.mu_kernel = create_kernel(mu_assignments, **kernel_parameters).compile()
+
+ force_rhs = force_from_phi_and_mu(self.phi_field.center_vector, dim=dh.dim, mu=self.mu_field.center_vector)
+ force_rhs = [discretize_spatial(e, dx=1, stencil=fd_stencils_forth_order_isotropic) for e in force_rhs]
+ force_assignments = [Assignment(lhs, rhs)
+ for lhs, rhs in zip(self.force_field.center_vector, force_rhs)]
+ self.force_kernel = create_kernel(force_assignments, **kernel_parameters).compile()
+
+ self.ch_lb_kernels = []
+ for i, (src, dst) in enumerate(self.ch_pdfs):
+ ch_method = cahn_hilliard_lb_method(stencil, self.mu_field(i),
+
+ relaxation_rate=cahn_hilliard_relaxation_rates[i],
+ gamma=cahn_hilliard_gammas[i])
+ opt = optimization.copy()
+ opt['symbolic_field'] = src
+ opt['symbolic_temporary_field'] = dst
+ kernel = create_lb_function(lb_method=ch_method, optimization=opt,
+ velocity_input=self.vel_field.center_vector,
+ output={'density': self.phi_field(i)})
+ self.ch_lb_kernels.append(kernel)
+
+ opt = optimization.copy()
+ opt['symbolic_field'] = self.hydro_pdfs[0]
+ opt['symbolic_temporary_field'] = self.hydro_pdfs[1]
+ self.hydro_lb_kernel = create_lb_function(stencil=stencil, relaxation_rate=hydro_dynamic_relaxation_rate,
+ force=self.force_field.center_vector,
+ output={'velocity': self.vel_field}, optimization=opt)
+
+ # Setter Kernels
+ self.init_kernels = []
+ for i in range(len(order_parameters)):
+ ch_method = self.ch_lb_kernels[i].method
+ init_assign = pdf_initialization_assignments(lb_method=ch_method,
+ density=self.phi_field.center_vector[i],
+ velocity=self.vel_field.center_vector,
+ pdfs=self.ch_pdfs[i][0].center_vector)
+ init_kernel = create_kernel(init_assign, **kernel_parameters).compile()
+ self.init_kernels.append(init_kernel)
+
+ init_assign = pdf_initialization_assignments(lb_method=self.hydro_lb_kernel.method, density=1,
+ velocity=self.vel_field.center_vector,
+ pdfs=self.hydro_pdfs[0].center_vector)
+ self.init_kernels.append(create_kernel(init_assign, **kernel_parameters).compile())
+
+ # Sync functions
+ self.phi_sync = dh.synchronization_function([self.phi_field.name])
+ self.mu_sync = dh.synchronization_function([self.mu_field.name])
+ self.pdf_sync = dh.synchronization_function([self.hydro_pdfs[0].name] +
+ [src.name for src, _ in self.ch_pdfs])
+
+ self.reset()
+
+ def reset(self):
+ dh = self.data_handling
+ dh.fill(self.vel_field.name, 0)
+ dh.fill(self.mu_field.name, 0)
+ dh.fill(self.force_field.name, 0)
+
+ def set_pdf_fields_from_macroscopic_values(self):
+ self.reset()
+ for k in self.init_kernels:
+ self.data_handling.run_kernel(k)
+
+ def set_concentration(self, slice_obj, concentration):
+ phi = np.array(concentration)
+
+ for b in self.data_handling.iterate(slice_obj):
+ for i in range(phi.shape[-1]):
+ b[self.phi_field.name][..., i] = phi[i]
+
+ def set_single_concentration(self, slice_obj, phase_idx, value=1):
+ num_phases = self.phi_field.index_shape[0]
+ zeros = [0] * num_phases
+ zeros[phase_idx] = value
+ self.set_concentration(slice_obj, zeros)
+
+ def smooth(self, sigma=2):
+ from scipy.ndimage.filters import gaussian_filter
+ dh = self.data_handling
+
+ for block in dh.iterate(ghost_layers=True):
+ c_arr = block[self.phi_field.name]
+ for i in range(self.phi_field.index_shape[0]):
+ gaussian_filter(c_arr[..., i], sigma=sigma, output=c_arr[..., i])
+
+ def time_step(self):
+ dh = self.data_handling
+
+ self.phi_sync()
+ dh.run_kernel(self.mu_kernel)
+
+ self.mu_sync()
+ dh.run_kernel(self.force_kernel)
+
+ self.pdf_sync()
+ dh.run_kernel(self.hydro_lb_kernel)
+ dh.swap(self.hydro_pdfs[0].name, self.hydro_pdfs[1].name)
+
+ for ch_kernel, (src, dst) in zip(self.ch_lb_kernels, self.ch_pdfs):
+ dh.run_kernel(ch_kernel)
+ dh.swap(src.name, dst.name)
--
GitLab