diff --git a/CHT/flow-over-plate/buoyantPimpleFoam-fenics/Solid/heat.py b/CHT/flow-over-plate/buoyantPimpleFoam-fenics/Solid/heat.py index 244efffbf..b2cce5ff1 100644 --- a/CHT/flow-over-plate/buoyantPimpleFoam-fenics/Solid/heat.py +++ b/CHT/flow-over-plate/buoyantPimpleFoam-fenics/Solid/heat.py @@ -6,7 +6,7 @@ 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 + near, MeshFunction, MPI from fenicsprecice import Adapter import numpy as np @@ -140,8 +140,17 @@ def determine_heat_flux(V_g, u, k, flux): t = 0 u_D.t = t + dt +# mark mesh w.r.t ranks +ranks = File("Solid/VTK/ranks%s.pvd.pvd" % precice.get_participant_name()) +mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim()) +mesh_rank.set_all(MPI.rank(MPI.comm_world)) +mesh_rank.rename("myRank", "") +ranks << mesh_rank + +# Create output file file_out = File("Solid/VTK/%s.pvd" % precice.get_participant_name()) file_out << u_n + print("output vtk for time = {}".format(float(t))) n = 0 diff --git a/CHT/flow-over-plate/buoyantPimpleFoam-fenics/precice-config.xml b/CHT/flow-over-plate/buoyantPimpleFoam-fenics/precice-config.xml index 87df65b75..df507989f 100644 --- a/CHT/flow-over-plate/buoyantPimpleFoam-fenics/precice-config.xml +++ b/CHT/flow-over-plate/buoyantPimpleFoam-fenics/precice-config.xml @@ -6,7 +6,7 @@ - + @@ -23,8 +23,10 @@ + + @@ -32,8 +34,7 @@ - - + @@ -41,14 +42,13 @@ - - + + - - + diff --git a/FSI/cylinderFlap/OpenFOAM-FEniCS/precice-config.xml b/FSI/cylinderFlap/OpenFOAM-FEniCS/precice-config.xml index 5b333dede..5c5b3ca2d 100644 --- a/FSI/cylinderFlap/OpenFOAM-FEniCS/precice-config.xml +++ b/FSI/cylinderFlap/OpenFOAM-FEniCS/precice-config.xml @@ -6,7 +6,7 @@ - + @@ -30,15 +30,15 @@ - - + + - + diff --git a/FSI/flap_perp/OpenFOAM-FEniCS/Solid/perp-flap.py b/FSI/flap_perp/OpenFOAM-FEniCS/Solid/perp-flap.py index 41dfeff4b..4b1891d17 100644 --- a/FSI/flap_perp/OpenFOAM-FEniCS/Solid/perp-flap.py +++ b/FSI/flap_perp/OpenFOAM-FEniCS/Solid/perp-flap.py @@ -1,7 +1,7 @@ # Import required libs from fenics import Constant, Function, AutoSubDomain, RectangleMesh, VectorFunctionSpace, interpolate, \ - TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \ - Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system + TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, SubDomain, \ + Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, PointSource, assemble_system, MPI, MeshFunction import dolfin from ufl import nabla_div @@ -11,17 +11,22 @@ from enum import Enum -# define the two kinds of boundary: clamped and coupling Neumann Boundary -def clamped_boundary(x, on_boundary): - return on_boundary and abs(x[1]) < tol +class clampedBoundary(SubDomain): + def inside(self, x, on_boundary): + tol = 1E-14 + if on_boundary and abs(x[1]) < tol: + return True + else: + return False -def neumann_boundary(x, on_boundary): - """ - determines whether a node is on the coupling boundary - - """ - return on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol) +class neumannBoundary(SubDomain): + def inside(self, x, on_boundary): + tol = 1E-14 + if on_boundary and ((abs(x[1] - 1) < tol) or abs(abs(x[0]) - W / 2) < tol): + return True + else: + return False # Geometry and material properties @@ -46,9 +51,6 @@ def neumann_boundary(x, on_boundary): # create Function Space V = VectorFunctionSpace(mesh, 'P', 2) -# BCs -tol = 1E-14 - # Trial and Test Functions du = TrialFunction(V) v = TestFunction(V) @@ -64,24 +66,19 @@ def neumann_boundary(x, on_boundary): f_N_function = interpolate(Expression(("1", "0"), degree=1), V) u_function = interpolate(Expression(("0", "0"), degree=1), V) -coupling_boundary = AutoSubDomain(neumann_boundary) - precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json") -clamped_boundary_domain = AutoSubDomain(clamped_boundary) -force_boundary = AutoSubDomain(neumann_boundary) - # Initialize the coupling interface # Function space V is passed twice as both read and write functions are defined using the same space -precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, - fixed_boundary=clamped_boundary_domain) +precice_dt = precice.initialize(neumannBoundary(), read_function_space=V, write_object=V, + fixed_boundary=clampedBoundary()) fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied # fenics_dt = 0.02 # if fenics_dt < precice_dt, subcycling is applied dt = Constant(np.min([precice_dt, fenics_dt])) # clamp the beam at the bottom -bc = DirichletBC(V, Constant((0, 0)), clamped_boundary) +bc = DirichletBC(V, Constant((0, 0)), clampedBoundary()) # alpha method parameters alpha_m = Constant(0.2) @@ -179,11 +176,18 @@ def avg(x_old, x_new, alpha): u_tip.append(0.0) E_ext = 0 +# mark mesh w.r.t ranks +mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim()) +mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0) +mesh_rank.rename("myRank", "") + displacement_out = File("Solid/FSI-S/u_fsi.pvd") +ranks = File("Solid/FSI-S/ranks%s.pvd" % precice.get_participant_name()) u_n.rename("Displacement", "") u_np1.rename("Displacement", "") displacement_out << u_n +ranks << mesh_rank while precice.is_coupling_ongoing(): diff --git a/FSI/flap_perp/OpenFOAM-FEniCS/precice-config.xml b/FSI/flap_perp/OpenFOAM-FEniCS/precice-config.xml index 385a6288b..3f42e25ce 100644 --- a/FSI/flap_perp/OpenFOAM-FEniCS/precice-config.xml +++ b/FSI/flap_perp/OpenFOAM-FEniCS/precice-config.xml @@ -6,7 +6,7 @@ - + @@ -30,15 +30,15 @@ - - + + - + diff --git a/HT/partitioned-heat/fenics-fenics/heat.py b/HT/partitioned-heat/fenics-fenics/heat.py index 8b98af863..8a0a4757d 100644 --- a/HT/partitioned-heat/fenics-fenics/heat.py +++ b/HT/partitioned-heat/fenics-fenics/heat.py @@ -26,7 +26,7 @@ from __future__ import print_function, division from fenics import Function, FunctionSpace, Expression, Constant, DirichletBC, TrialFunction, TestFunction, \ - File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace + File, solve, lhs, rhs, grad, inner, dot, dx, ds, interpolate, VectorFunctionSpace, MeshFunction, MPI from fenicsprecice import Adapter from errorcomputation import compute_errors from my_enums import ProblemType, Subcycling @@ -174,16 +174,26 @@ def determine_gradient(V_g, u, flux): u_ref = interpolate(u_D, V) u_ref.rename("reference", " ") +# mark mesh w.r.t ranks +mesh_rank = MeshFunction("size_t", mesh, mesh.topology().dim()) +if problem is ProblemType.NEUMANN: + mesh_rank.set_all(MPI.rank(MPI.comm_world) + 4) +else: + mesh_rank.set_all(MPI.rank(MPI.comm_world) + 0) +mesh_rank.rename("myRank", "") + # Generating output files temperature_out = File("out/%s.pvd" % precice.get_participant_name()) ref_out = File("out/ref%s.pvd" % precice.get_participant_name()) error_out = File("out/error%s.pvd" % precice.get_participant_name()) +ranks = File("out/ranks%s.pvd" % precice.get_participant_name()) # output solution and reference solution at t=0, n=0 n = 0 print('output u^%d and u_ref^%d' % (n, n)) temperature_out << u_n ref_out << u_ref +ranks << mesh_rank error_total, error_pointwise = compute_errors(u_n, u_ref, V) error_out << error_pointwise diff --git a/HT/partitioned-heat/fenics-fenics/precice-config.xml b/HT/partitioned-heat/fenics-fenics/precice-config.xml index 5c33d9632..129a7fcf0 100644 --- a/HT/partitioned-heat/fenics-fenics/precice-config.xml +++ b/HT/partitioned-heat/fenics-fenics/precice-config.xml @@ -26,14 +26,15 @@ - - + + + @@ -43,9 +44,9 @@ - + - +