Use same method for flux computation like in HT example. (#114)

Author: BenjaminRodenberg

CHT/flow-over-plate/buoyantPimpleFoam-fenics/Solid/heat.py

@@ -3,8 +3,9 @@
 """
 
 from __future__ import print_function, division
-from fenics import Function, SubDomain, RectangleMesh, BoxMesh, FunctionSpace, Point, Expression, Constant, DirichletBC, \
-    TrialFunction, TestFunction, File, solve, plot, lhs, rhs, grad, inner, dot, dx, ds, assemble, interpolate, project, \
+from fenics import Function, SubDomain, RectangleMesh, BoxMesh, FunctionSpace, VectorFunctionSpace, Point, \
+    Expression, Constant, DirichletBC, \
+    TrialFunction, TestFunction, File, solve, plot, lhs, rhs, grad, inner, dot, dx, ds, interpolate, project, \
     near
 from fenicsadapter import Adapter
 import numpy as np
@@ -59,27 +60,21 @@ def inside(self, x, on_boundary):
             return False
 
 
-def fluxes_from_temperature_full_domain(f, v_vec, k):
-    """Computes flux from weak form (see p.3 in Toselli, Andrea, and Olof
-    Widlund. Domain decomposition methods-algorithms and theory. Vol. 34.
-    Springer Science & Business Media, 2006.).
-
-    :param f: weak form with known u^{n+1}
-    :param v_vec: vector function space
+def determine_heat_flux(V_g, u, k, flux):
+    """
+    compute flux following http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html#tut-poisson-gradu
+    :param V_g: Vector function space
+    :param u: solution where gradient is to be determined
     :param k: thermal conductivity
-    :return: fluxes function
+    :param flux: returns calculated flux into this value
     """
-    fluxes_vector = assemble(f)  # assemble weak form -> evaluate integral
-    v = TestFunction(v_vec)
-    fluxes = Function(v_vec)  # create function for flux
-    area = assemble(v * ds).get_local()
-    for i in range(area.shape[0]):
-        if area[i] != 0:  # put weight from assemble on function
-            fluxes.vector()[i] = - k * fluxes_vector[i] / area[i]  # scale by surface area
-        else:
-            assert (abs(fluxes_vector[i]) < 1E-9)  # for non surface parts, we expect zero flux
-            fluxes.vector()[i] = - k * fluxes_vector[i]
-    return fluxes
+
+    w = TrialFunction(V_g)
+    v = TestFunction(V_g)
+
+    a = inner(w, v) * dx
+    L = inner(-k * grad(u), v) * dx
+    solve(a == L, flux)
 
 
 # Create mesh and define function space
@@ -99,14 +94,15 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 
 mesh = RectangleMesh(p0, p1, nx, ny)
 V = FunctionSpace(mesh, 'P', 1)
+V_g = VectorFunctionSpace(mesh, 'P', 1)
 
 alpha = 1  # m^2/s, https://en.wikipedia.org/wiki/Thermal_diffusivity
 k = 100  # kg * m / s^3 / K, https://en.wikipedia.org/wiki/Thermal_conductivity
 
 # Define boundary condition
 u_D = Constant('310')
 u_D_function = interpolate(u_D, V)
-# Define flux in x direction on coupling interface (grad(u_D) in normal direction)
+# Define flux in y direction on coupling interface (grad(u_D) in normal direction)
 f_N = Constant('0')
 f_N_function = interpolate(f_N, V)
 
@@ -120,7 +116,7 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 # Adapter definition and initialization
 precice = Adapter(adapter_config_filename="precice-adapter-config.json")
 
-precice_dt = precice.initialize(coupling_boundary, mesh, V)
+precice_dt = precice.initialize(coupling_boundary, mesh, V, write_function=f_N_function)
 
 # Create a FEniCS Expression to define and control the coupling boundary values
 coupling_expression = precice.create_coupling_expression()
@@ -142,7 +138,6 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 
 # Time-stepping
 u_np1 = Function(V)
-F_known_u = u_np1 * v / dt * dx + alpha * dot(grad(u_np1), grad(v)) * dx - u_n * v / dt * dx
 t = 0
 u_D.t = t + dt
 
@@ -151,6 +146,9 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
 print("output vtk for time = {}".format(float(t)))
 n = 0
 
+fluxes = Function(V_g)
+fluxes.rename("Fluxes", "")
+
 while precice.is_coupling_ongoing():
 
     if precice.is_action_required(precice.action_write_iteration_checkpoint()):  # write checkpoint
@@ -167,8 +165,9 @@ def fluxes_from_temperature_full_domain(f, v_vec, k):
     solve(a == L, u_np1, bcs)
 
     # Dirichlet problem obtains flux from solution and sends flux on boundary to Neumann problem
-    fluxes = fluxes_from_temperature_full_domain(F_known_u, V, k)
-    precice.write_data(fluxes)
+    determine_heat_flux(V_g, u_np1, k, fluxes)
+    fluxes_y = fluxes.sub(1)  # only exchange y component of flux.
+    precice.write_data(fluxes_y)
 
     precice_dt = precice.advance(dt(0))