From 4af616ff56e5ec2e58536543e15d45144d4ffc39 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Wed, 5 Dec 2018 12:15:02 +0100
Subject: [PATCH] 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
---
phasefield/analytical.py | 22 ++--------------------
phasefield/kerneleqs.py | 4 ++--
phasefield/phasefieldstep.py | 19 ++++---------------
3 files changed, 8 insertions(+), 37 deletions(-)
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 04141d74..93a6234a 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -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
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index fce9bbed..75f421d9 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -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)
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index 706e0c07..31b070d9 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -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),
--
GitLab