Formatting file paritioned-heat-conduction/nutils/heat.py

Author: IshaanDesai

partitioned-heat-conduction/nutils/heat.py

@@ -19,7 +19,8 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
 
     # define the Nutils mesh
     domain, geom = mesh.rectilinear([x_grid, y_grid])
-    coupling_boundary = domain.boundary['right' if side == 'Dirichlet' else 'left']
+    coupling_boundary = domain.boundary['right' if side ==
+                                        'Dirichlet' else 'left']
     coupling_sample = coupling_boundary.sample('gauss', degree=degree * 2)
 
     # Nutils namespace
@@ -37,12 +38,14 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
     ns.readfunc = 'readbasis_n ?readdata_n'
 
     # define the weak form
-    res = domain.integral('(basis_n dudt - basis_n f + basis_n,i u_,i) d:x' @ ns, degree=degree * 2)
+    res = domain.integral(
+        '(basis_n dudt - basis_n f + basis_n,i u_,i) d:x' @ ns,
+        degree=degree * 2)
 
     # set boundary conditions at non-coupling boundaries
     # top and bottom boundary are non-coupling for both sides
-    sqr = domain.boundary['top,bottom,left' if side == 'Dirichlet'
-                     else 'top,bottom,right'].integral('(u - uexact)^2 d:x' @ ns, degree=degree * 2)
+    sqr = domain.boundary['top,bottom,left' if side == 'Dirichlet' else 'top,bottom,right'].integral(
+        '(u - uexact)^2 d:x' @ ns, degree=degree * 2)
 
     if side == 'Dirichlet':
         sqr += coupling_sample.integral('(u - readfunc)^2 d:x' @ ns)
@@ -52,11 +55,20 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
     # preCICE setup
     interface = precice.Interface(side, "../precice-config.xml", 0, 1)
     mesh_id = interface.get_mesh_id(side + "-Mesh")
-    vertex_ids = interface.set_mesh_vertices(mesh_id, coupling_sample.eval(ns.x))
-    precice_write = functools.partial(interface.write_block_scalar_data,
-        interface.get_data_id("Temperature" if side == "Neumann" else "Heat-Flux", mesh_id), vertex_ids)
-    precice_read = functools.partial(interface.read_block_scalar_data,
-        interface.get_data_id("Heat-Flux" if side == "Neumann" else "Temperature", mesh_id), vertex_ids)
+    vertex_ids = interface.set_mesh_vertices(
+        mesh_id, coupling_sample.eval(ns.x))
+    precice_write = functools.partial(
+        interface.write_block_scalar_data,
+        interface.get_data_id(
+            "Temperature" if side == "Neumann" else "Heat-Flux",
+            mesh_id),
+        vertex_ids)
+    precice_read = functools.partial(
+        interface.read_block_scalar_data,
+        interface.get_data_id(
+            "Heat-Flux" if side == "Neumann" else "Temperature",
+            mesh_id),
+        vertex_ids)
 
     # helper functions to project heat flux to coupling boundary
     if side == 'Dirichlet':
@@ -70,11 +82,18 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
         # While the latter still contains the problematic unbounded term, we
         # can use the fact that the flux is a known value at the top and bottom
         # via the Dirichlet boundary condition, and impose it as constraints.
-        rightsqr = domain.boundary['right'].integral('flux^2 d:x' @ ns, degree=degree*2)
+        rightsqr = domain.boundary['right'].integral(
+            'flux^2 d:x' @ ns, degree=degree * 2)
         rightcons = solver.optimize('fluxdofs', rightsqr, droptol=1e-10)
         # rightcons is NaN in dofs that are NOT supported on the right boundary
-        fluxsqr = domain.boundary['right'].boundary['top,bottom'].integral('(flux - uexact_,0)^2 d:x' @ ns, degree=degree*2)
-        fluxcons = solver.optimize('fluxdofs', fluxsqr, droptol=1e-10, constrain=np.choose(np.isnan(rightcons), [np.nan, 0.]))
+        fluxsqr = domain.boundary['right'].boundary['top,bottom'].integral(
+            '(flux - uexact_,0)^2 d:x' @ ns, degree=degree * 2)
+        fluxcons = solver.optimize('fluxdofs',
+                                   fluxsqr,
+                                   droptol=1e-10,
+                                   constrain=np.choose(np.isnan(rightcons),
+                                                       [np.nan,
+                                                        0.]))
         # fluxcons is NaN in dofs that are supported on ONLY the right boundary
         fluxres = coupling_sample.integral('basis_n flux d:x' @ ns) - res
 
@@ -100,7 +119,14 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
         # generate output
         x, u, uexact = bezier.eval(['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t)
         with treelog.add(treelog.DataLog()):
-            export.vtk(side + "-" + str(istep), bezier.tri, x, Temperature=u, reference=uexact)
+            export.vtk(
+                side +
+                "-" +
+                str(istep),
+                bezier.tri,
+                x,
+                Temperature=u,
+                reference=uexact)
 
         if not interface.is_coupling_ongoing():
             break
@@ -110,9 +136,11 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
             readdata = precice_read()
 
         # save checkpoint
-        if interface.is_action_required(precice.action_write_iteration_checkpoint()):
+        if interface.is_action_required(
+                precice.action_write_iteration_checkpoint()):
             checkpoint = lhs, t, istep
-            interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())
+            interface.mark_action_fulfilled(
+                precice.action_write_iteration_checkpoint())
 
         # prepare next timestep
         lhs0 = lhs
@@ -121,16 +149,27 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
         t += dt
 
         # update (time-dependent) boundary condition
-        cons = solver.optimize('lhs', sqr, droptol=1e-15, arguments=dict(t=t, readdata=readdata))
+        cons = solver.optimize(
+            'lhs',
+            sqr,
+            droptol=1e-15,
+            arguments=dict(
+                t=t,
+                readdata=readdata))
 
         # solve nutils timestep
-        lhs = solver.solve_linear('lhs', res, constrain=cons, arguments=dict(lhs0=lhs0, dt=dt, t=t, readdata=readdata))
+        lhs = solver.solve_linear(
+            'lhs', res, constrain=cons, arguments=dict(
+                lhs0=lhs0, dt=dt, t=t, readdata=readdata))
 
         # write data to interface
         if interface.is_write_data_required(dt):
             if side == 'Dirichlet':
-                fluxdofs = solver.solve_linear('fluxdofs', fluxres, arguments=dict(lhs0=lhs0, lhs=lhs, dt=dt, t=t), constrain=fluxcons)
-                write_data = coupling_sample.eval('flux' @ ns, fluxdofs=fluxdofs)
+                fluxdofs = solver.solve_linear(
+                    'fluxdofs', fluxres, arguments=dict(
+                        lhs0=lhs0, lhs=lhs, dt=dt, t=t), constrain=fluxcons)
+                write_data = coupling_sample.eval(
+                    'flux' @ ns, fluxdofs=fluxdofs)
             else:
                 write_data = coupling_sample.eval('u' @ ns, lhs=lhs)
             precice_write(write_data)
@@ -139,9 +178,11 @@ def main(side='Dirichlet', n=10, degree=1, timestep=.1, alpha=3., beta=1.3):
         precice_dt = interface.advance(dt)
 
         # read checkpoint if required
-        if interface.is_action_required(precice.action_read_iteration_checkpoint()):
+        if interface.is_action_required(
+                precice.action_read_iteration_checkpoint()):
             lhs, t, istep = checkpoint
-            interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())
+            interface.mark_action_fulfilled(
+                precice.action_read_iteration_checkpoint())
 
     interface.finalize()