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

Removed transformation before pressure tensor computation

- conceptionally not required: force needs to be in transformed coords
- for 3 phase model it makes no difference after bug was fixed in
  extract_gammas
parent 3a78346c
Branches
Tags
No related merge requests found
......@@ -343,22 +343,9 @@ def pressure_tensor_interface_component_new(free_energy, order_parameters, dim,
return result
def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix=None,
def pressure_tensor_from_free_energy(free_energy, order_parameters, dim,
include_bulk=True, include_interface=True):
def transform(f, matrix, from_symbols, to_symbols):
transformed_vals = matrix * sp.Matrix(to_symbols)
substitutions = {a: b for a, b in zip(from_symbols, transformed_vals)}
return f.subs(substitutions)
c = sp.symbols("C_:{}".format(len(order_parameters)))
if transformation_matrix:
phi_to_c = partial(transform, matrix=transformation_matrix, from_symbols=order_parameters, to_symbols=c)
free_energy = expand_diff_full(phi_to_c(free_energy))
op = c
else:
op = order_parameters
op = order_parameters
def get_entry(i, j):
p_if = pressure_tensor_interface_component(free_energy, op, dim, i, j) if include_interface else 0
......@@ -369,11 +356,6 @@ def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transfo
return sp.expand(p_if + p_b)
result = sp.Matrix(dim, dim, get_entry)
if transformation_matrix:
c_to_phi = partial(transform, matrix=transformation_matrix.inv(), from_symbols=c, to_symbols=order_parameters)
result = result.applyfunc(lambda e: expand_diff_full(c_to_phi(e)))
return result
......
......@@ -29,10 +29,10 @@ def force_kernel_using_mu(force_field, phi_field, mu_field, dx=1, discretization
def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_tensor_field,
transformation_matrix=None, dx=1, discretization='standard'):
dx=1, discretization='standard'):
dim = phi_field.spatial_dimensions
p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix)
p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim)
p = p.subs({op: phi_field(i) for i, op in enumerate(order_parameters)})
index_map = symmetric_tensor_linearization(dim)
......
......@@ -6,7 +6,7 @@ import numpy as np
from lbmpy.lbstep import LatticeBoltzmannStep
from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
from lbmpy.phasefield.kerneleqs import mu_kernel, CahnHilliardFDStep, pressure_tensor_kernel, \
force_kernel_using_pressure_tensor, pressure_tensor_kernel_pbs
force_kernel_using_pressure_tensor
from pystencils import create_kernel, create_data_handling
from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, symmetric_tensor_linearization
from pystencils.boundaries.boundaryhandling import FlagInterface
......@@ -98,20 +98,9 @@ class PhaseFieldStep:
self.phi_sync = data_handling.synchronization_function([self.phi_field_name], target=target)
self.mu_eqs = mu_kernel(F, phi, self.phi_field, self.mu_field, dx)
pbs = False
if pbs:
self.pbs_field_name = name + "_pbs"
self.pbs_field = dh.add_array(self.pbs_field_name, gpu=gpu)
self.pressure_tensor_eqs = pressure_tensor_kernel_pbs(self.free_energy, order_parameters, self.phi_field,
self.pressure_tensor_field, self.pbs_field, dx=dx,
density_field=None, discretization=discretization)
# TODO get current density! not last one
# TODO call post-run on hydro-lbm before computing pbs to store the latest density
else:
self.pressure_tensor_eqs = pressure_tensor_kernel(self.free_energy, order_parameters,
self.phi_field, self.pressure_tensor_field, dx=dx,
transformation_matrix=transformation_matrix,
discretization=discretization)
self.pressure_tensor_eqs = pressure_tensor_kernel(self.free_energy, order_parameters,
self.phi_field, self.pressure_tensor_field, dx=dx,
discretization=discretization)
mu_and_pressure_tensor_eqs = self.mu_eqs + self.pressure_tensor_eqs
mu_and_pressure_tensor_eqs = apply_neumann_boundaries(mu_and_pressure_tensor_eqs)
self.mu_and_pressure_tensor_kernel = create_kernel(sympy_cse_on_assignment_list(mu_and_pressure_tensor_eqs),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment