Skip to content
Snippets Groups Projects
Commit 988e05a4 authored by Martin Bauer's avatar Martin Bauer
Browse files

lbm phase field notebooks refactoring

parent 72e4fc40
No related branches found
No related tags found
No related merge requests found
...@@ -170,26 +170,28 @@ def separate_into_bulk_and_interface(free_energy): ...@@ -170,26 +170,28 @@ def separate_into_bulk_and_interface(free_energy):
def analytic_interface_profile(x, interface_width=interface_width_symbol): def analytic_interface_profile(x, interface_width=interface_width_symbol):
"""Analytic expression for a 1D interface normal to x with given interface width """Analytic expression for a 1D interface normal to x with given interface width.
The following doctest shows that the returned analytical solution is indeed a solution of the ODE that we The following doctest shows that the returned analytical solution is indeed a solution of the ODE that we
get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order
parameter, i.e. at a transition between two phases. parameter, i.e. at a transition between two phases.
>>> numPhases = 4
>>> x, phi = sp.Symbol("x"), symbolic_order_parameters(numPhases-1) Examples:
>>> F = free_energy_functional_n_phases(order_parameters=phi) >>> numPhases = 4
>>> mu = chemical_potentials_from_free_energy(F) >>> x, phi = sp.Symbol("x"), symbolic_order_parameters(numPhases-1)
>>> mu0 = mu[0].subs({p: 0 for p in phi[1:]}) # mu[0] as function of one order parameter only >>> F = free_energy_functional_n_phases(order_parameters=phi)
>>> solution = analytic_interface_profile(x) >>> mu = chemical_potentials_from_free_energy(F)
>>> solutionSubstitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) } >>> mu0 = mu[0].subs({p: 0 for p in phi[1:]}) # mu[0] as function of one order parameter only
>>> sp.expand(mu0.subs(solutionSubstitution)) # inserting solution should solve the mu_0=0 equation >>> solution = analytic_interface_profile(x)
0 >>> solutionSubstitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) }
>>> sp.expand(mu0.subs(solutionSubstitution)) # inserting solution should solve the mu_0=0 equation
0
""" """
return (1 + sp.tanh(x / (2 * interface_width))) / 2 return (1 + sp.tanh(x / (2 * interface_width))) / 2
def chemical_potentials_from_free_energy(free_energy, order_parameters=None): def chemical_potentials_from_free_energy(free_energy, order_parameters=None):
"""Computes chemical potentials as functional derivative of free energy""" """Computes chemical potentials as functional derivative of free energy."""
symbols = free_energy.atoms(sp.Symbol) symbols = free_energy.atoms(sp.Symbol)
if order_parameters is None: if order_parameters is None:
order_parameters = [s for s in symbols if s.name.startswith(orderParameterSymbolName)] order_parameters = [s for s in symbols if s.name.startswith(orderParameterSymbolName)]
......
...@@ -13,26 +13,26 @@ def plot_status(phase_field_step, from_x=None, to_x=None): ...@@ -13,26 +13,26 @@ def plot_status(phase_field_step, from_x=None, to_x=None):
dh = phase_field_step.data_handling dh = phase_field_step.data_handling
num_phases = phase_field_step.numOrderParameters num_phases = phase_field_step.num_order_parameters
plt.subplot(1, 3, 1) plt.subplot(1, 3, 1)
plt.title('φ') plt.title('φ')
phi_name = phase_field_step.phiFieldName phi_name = phase_field_step.phi_field_name
for i in range(num_phases): for i in range(num_phases):
plt.plot(dh.gather_array(phi_name, make_slice[from_x:to_x, 0, i]), marker='x', label='φ_%d' % (i,)) plt.plot(dh.gather_array(phi_name, make_slice[from_x:to_x, 0, i]), marker='x', label='φ_%d' % (i,))
plt.legend() plt.legend()
plt.subplot(1, 3, 2) plt.subplot(1, 3, 2)
plt.title("μ") plt.title("μ")
mu_name = phase_field_step.muFieldName mu_name = phase_field_step.mu_field_name
for i in range(num_phases): for i in range(num_phases):
plt.plot(dh.gather_array(mu_name, make_slice[from_x:to_x, 0, i]), marker='x', label='μ_%d' % (i,)) plt.plot(dh.gather_array(mu_name, make_slice[from_x:to_x, 0, i]), marker='x', label='μ_%d' % (i,))
plt.legend() plt.legend()
plt.subplot(1, 3, 3) plt.subplot(1, 3, 3)
plt.title("Force and Velocity") plt.title("Force and Velocity")
plt.plot(dh.gather_array(phase_field_step.forceFieldName, make_slice[from_x:to_x, 0, 0]), label='F', marker='x') plt.plot(dh.gather_array(phase_field_step.force_field_name, make_slice[from_x:to_x, 0, 0]), label='F', marker='x')
plt.plot(dh.gather_array(phase_field_step.velFieldName, make_slice[from_x:to_x, 0, 0]), label='u', marker='v') plt.plot(dh.gather_array(phase_field_step.vel_field_name, make_slice[from_x:to_x, 0, 0]), label='u', marker='v')
plt.legend() plt.legend()
...@@ -59,14 +59,14 @@ def init_sharp_interface(pf_step, phase_idx, x1=None, x2=None, inverse=False): ...@@ -59,14 +59,14 @@ def init_sharp_interface(pf_step, phase_idx, x1=None, x2=None, inverse=False):
if x2 is None: if x2 is None:
x2 = 3 * x1 x2 = 3 * x1
if phase_idx >= pf_step.numOrderParameters: if phase_idx >= pf_step.num_order_parameters:
return return
for b in pf_step.data_handling.iterate(): for b in pf_step.data_handling.iterate():
x = b.cell_index_arrays[0] x = b.cell_index_arrays[0]
mid = np.logical_and(x1 < x, x < x2) mid = np.logical_and(x1 < x, x < x2)
phi = b[pf_step.phiFieldName] phi = b[pf_step.phi_field_name]
val1, val2 = (1, 0) if inverse else (0, 1) val1, val2 = (1, 0) if inverse else (0, 1)
phi[..., phase_idx].fill(val1) phi[..., phase_idx].fill(val1)
...@@ -94,7 +94,7 @@ def tanh_test(pf_step, phase0, phase1, expected_interface_width=1, time_steps=10 ...@@ -94,7 +94,7 @@ def tanh_test(pf_step, phase0, phase1, expected_interface_width=1, time_steps=10
domain_size = pf_step.data_handling.shape domain_size = pf_step.data_handling.shape
pf_step.reset() pf_step.reset()
pf_step.data_handling.fill(pf_step.phiFieldName, 0) pf_step.data_handling.fill(pf_step.phi_field_name, 0)
init_sharp_interface(pf_step, phase_idx=phase0, inverse=False) init_sharp_interface(pf_step, phase_idx=phase0, inverse=False)
init_sharp_interface(pf_step, phase_idx=phase1, inverse=True) init_sharp_interface(pf_step, phase_idx=phase1, inverse=True)
pf_step.set_pdf_fields_from_macroscopic_values() pf_step.set_pdf_fields_from_macroscopic_values()
...@@ -140,14 +140,14 @@ def galilean_invariance_test(pf_step, velocity=0.05, rounds=3, phase0=0, phase1= ...@@ -140,14 +140,14 @@ def galilean_invariance_test(pf_step, velocity=0.05, rounds=3, phase0=0, phase1=
print("Velocity:", velocity, " Time steps for round:", round_time_steps) print("Velocity:", velocity, " Time steps for round:", round_time_steps)
pf_step.reset() pf_step.reset()
pf_step.data_handling.fill(pf_step.phiFieldName, 0) pf_step.data_handling.fill(pf_step.phi_field_name, 0)
init_sharp_interface(pf_step, phase_idx=phase0, inverse=False) init_sharp_interface(pf_step, phase_idx=phase0, inverse=False)
init_sharp_interface(pf_step, phase_idx=phase1, inverse=True) init_sharp_interface(pf_step, phase_idx=phase1, inverse=True)
pf_step.set_pdf_fields_from_macroscopic_values() pf_step.set_pdf_fields_from_macroscopic_values()
print("Running", init_time_steps, "initial time steps") print("Running", init_time_steps, "initial time steps")
pf_step.run(init_time_steps) pf_step.run(init_time_steps)
pf_step.data_handling.fill(pf_step.velFieldName, velocity, value_idx=0) pf_step.data_handling.fill(pf_step.vel_field_name, velocity, value_idx=0)
pf_step.set_pdf_fields_from_macroscopic_values() pf_step.set_pdf_fields_from_macroscopic_values()
step_location = domain_size[0] // 4 step_location = domain_size[0] // 4
......
...@@ -30,11 +30,11 @@ class PhaseFieldStep: ...@@ -30,11 +30,11 @@ class PhaseFieldStep:
if data_handling is None: if data_handling is None:
data_handling = SerialDataHandling(domain_size, periodicity=True) data_handling = SerialDataHandling(domain_size, periodicity=True)
self.freeEnergy = free_energy self.free_energy = free_energy
self.concentrationToOrderParameter = concentration_to_order_parameters self.concentration_to_order_parameter = concentration_to_order_parameters
self.orderParametersToConcentrations = order_parameters_to_concentrations self.order_parameters_to_concentrations = order_parameters_to_concentrations
self.chemicalPotentials = chemical_potentials_from_free_energy(free_energy, order_parameters) self.chemical_potentials = chemical_potentials_from_free_energy(free_energy, order_parameters)
# ------------------ Adding arrays --------------------- # ------------------ Adding arrays ---------------------
gpu = target == 'gpu' gpu = target == 'gpu'
...@@ -61,7 +61,7 @@ class PhaseFieldStep: ...@@ -61,7 +61,7 @@ class PhaseFieldStep:
# ------------------ Creating kernels ------------------ # ------------------ Creating kernels ------------------
phi = tuple(self.phi_field(i) for i in range(len(order_parameters))) phi = tuple(self.phi_field(i) for i in range(len(order_parameters)))
F = self.freeEnergy.subs({old: new for old, new in zip(order_parameters, phi)}) F = self.free_energy.subs({old: new for old, new in zip(order_parameters, phi)})
if homogeneous_neumann_boundaries: if homogeneous_neumann_boundaries:
def apply_neumann_boundaries(eqs): def apply_neumann_boundaries(eqs):
...@@ -77,7 +77,7 @@ class PhaseFieldStep: ...@@ -77,7 +77,7 @@ class PhaseFieldStep:
# μ and pressure tensor update # μ and pressure tensor update
self.phiSync = data_handling.synchronization_function([self.phi_field_name], target=target) self.phiSync = data_handling.synchronization_function([self.phi_field_name], target=target)
self.muEqs = mu_kernel(F, phi, self.phi_field, self.mu_field, dx) self.muEqs = mu_kernel(F, phi, self.phi_field, self.mu_field, dx)
self.pressureTensorEqs = pressure_tensor_kernel(self.freeEnergy, order_parameters, self.pressureTensorEqs = pressure_tensor_kernel(self.free_energy, order_parameters,
self.phi_field, self.pressure_tensor_field, dx) self.phi_field, self.pressure_tensor_field, dx)
mu_and_pressure_tensor_eqs = self.muEqs + self.pressureTensorEqs mu_and_pressure_tensor_eqs = self.muEqs + self.pressureTensorEqs
mu_and_pressure_tensor_eqs = apply_neumann_boundaries(mu_and_pressure_tensor_eqs) mu_and_pressure_tensor_eqs = apply_neumann_boundaries(mu_and_pressure_tensor_eqs)
...@@ -225,8 +225,8 @@ class PhaseFieldStep: ...@@ -225,8 +225,8 @@ class PhaseFieldStep:
return self.hydroLbmStep.boundary_handling return self.hydroLbmStep.boundary_handling
def set_concentration(self, slice_obj, concentration): def set_concentration(self, slice_obj, concentration):
if self.concentrationToOrderParameter is not None: if self.concentration_to_order_parameter is not None:
phi = self.concentrationToOrderParameter(concentration) phi = self.concentration_to_order_parameter(concentration)
else: else:
phi = np.array(concentration) phi = np.array(concentration)
...@@ -255,7 +255,7 @@ class PhaseFieldStep: ...@@ -255,7 +255,7 @@ class PhaseFieldStep:
def concentration_slice(self, slice_obj=None): def concentration_slice(self, slice_obj=None):
phi = self.phi_slice(slice_obj) phi = self.phi_slice(slice_obj)
return phi if self.orderParametersToConcentrations is None else self.orderParametersToConcentrations(phi) return phi if self.order_parameters_to_concentrations is None else self.order_parameters_to_concentrations(phi)
def mu_slice(self, slice_obj=None): def mu_slice(self, slice_obj=None):
return self._get_slice(self.mu_field_name, slice_obj) return self._get_slice(self.mu_field_name, slice_obj)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment