Skip to content

Commit

Permalink
Update all FEniCS participants to be compatible with precice:develop …
Browse files Browse the repository at this point in the history
…(v3) (#333)

* Update all FEniCS participants to be compatible with precice:develop and fenics-adapter updates
* Change from thin plate splines mapping to general rbf mapping
* Update how dt is used.
* Store time in vtk files for better visualization.
* Partitioned heat conduction:
  * Fix config and case definition. (precice/precice#1610)
  * Use stricter error tolerance.
  • Loading branch information
IshaanDesai committed Sep 28, 2023
1 parent 0cd624e commit 63273dd
Show file tree
Hide file tree
Showing 13 changed files with 88 additions and 79 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,8 @@ def inside(self, x, on_boundary):
# Initialize preCICE
precice = fenicsprecice.Adapter(
adapter_config_filename="chemical-reaction-advection-diffusion.json")
precice_dt = precice.initialize(
coupling_subdomain=CouplingDomain(),
read_function_space=W)
precice.initialize(coupling_subdomain=CouplingDomain(), read_function_space=W)
precice_dt = precice.get_max_time_step_size()

flow_expr = precice.create_coupling_expression()

Expand Down Expand Up @@ -95,13 +94,14 @@ def inside(self, x, on_boundary):
# No implicit coupling
while precice.is_coupling_ongoing():

read_data = precice.read_data()
precice_dt = precice.get_max_time_step_size()
dt = np.min([default_dt, precice_dt])
read_data = precice.read_data(dt)
precice.update_coupling_expression(flow_expr, read_data)
flow_old.assign(flow)
flow.interpolate(flow_expr)
# If we add writing, do it here

dt = np.min([default_dt, precice_dt])
k.assign(1. / dt)

t += dt
Expand All @@ -128,6 +128,6 @@ def inside(self, x, on_boundary):
vtkfileC << u_C, t
vtkfileFlow << flow, t

precice_dt = precice.advance(dt)
precice.advance(dt)

precice.finalize()
9 changes: 4 additions & 5 deletions channel-transport-reaction/fluid-fenics/fluid.py
Original file line number Diff line number Diff line change
Expand Up @@ -118,10 +118,8 @@ def inside(self, x, on_boundary):


precice = fenicsprecice.Adapter(adapter_config_filename="fluid-config.json")
precice_dt = precice.initialize(
coupling_subdomain=CouplingDomain(),
write_object=u_)

precice.initialize(coupling_subdomain=CouplingDomain(), write_object=u_)
precice_dt = precice.get_max_time_step_size()

dt = np.min([default_dt, precice_dt])
k.assign(dt)
Expand Down Expand Up @@ -155,7 +153,8 @@ def inside(self, x, on_boundary):

vtkfile << u_

precice_dt = precice.advance(dt)
precice.advance(dt)
precice_dt = precice.get_max_time_step_size()
dt = np.min([default_dt, precice_dt])
k.assign(dt)

Expand Down
20 changes: 11 additions & 9 deletions elastic-tube-3d/solid-fenics/solid.py
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,8 @@ def neumann_boundary(x, on_boundary):
precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json")

# Initialize the coupling interface
precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)
precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary)
precice_dt = precice.get_max_time_step_size()

fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied
dt = Constant(np.min([precice_dt, fenics_dt]))
Expand Down Expand Up @@ -168,15 +169,18 @@ def avg(x_old, x_new, alpha):

u_n.rename("Displacement", "")
u_np1.rename("Displacement", "")
displacement_out << u_n
displacement_out << (u_n, t)

while precice.is_coupling_ongoing():

if precice.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint
if precice.requires_writing_checkpoint(): # write checkpoint
precice.store_checkpoint(u_n, t, n)

precice_dt = precice.get_max_time_step_size()
dt = Constant(np.min([precice_dt, fenics_dt]))

# read data from preCICE and get a new coupling expression
read_data = precice.read_data()
read_data = precice.read_data(dt)

# Update the point sources on the coupling boundary with the new read data
forces_x, forces_y, forces_z = precice.get_point_sources(read_data)
Expand All @@ -195,17 +199,15 @@ def avg(x_old, x_new, alpha):
assert (b is not b_forces)
solve(A, u_np1.vector(), b_forces)

dt = Constant(np.min([precice_dt, fenics_dt]))

# Write relative displacements to preCICE
u_delta.vector()[:] = u_np1.vector()[:] - u_n.vector()[:]
precice.write_data(u_delta)

# Call to advance coupling, also returns the optimum time step value
precice_dt = precice.advance(dt(0))
precice.advance(dt(0))

# Either revert to old step if timestep has not converged or move to next timestep
if precice.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint
if precice.requires_reading_checkpoint(): # roll back to checkpoint
u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
u_n.assign(u_cp)
t = t_cp
Expand All @@ -221,6 +223,6 @@ def avg(x_old, x_new, alpha):
displacement_out << (u_n, t)

# Plot tip displacement evolution
displacement_out << u_n
displacement_out << (u_n, t)

precice.finalize()
19 changes: 10 additions & 9 deletions flow-over-heated-plate/solid-fenics/solid.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,13 +115,14 @@ def determine_heat_flux(V_g, u, k, flux):
# Adapter definition and initialization
precice = Adapter(adapter_config_filename="precice-adapter-config.json")

precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V_flux_y)
precice.initialize(coupling_boundary, read_function_space=V, write_object=V_flux_y)

# Create a FEniCS Expression to define and control the coupling boundary values
coupling_expression = precice.create_coupling_expression()

# Assigning appropriate dt
dt = Constant(0)
precice_dt = precice.get_max_time_step_size()
dt.assign(np.min([fenics_dt, precice_dt]))

# Define variational problem
Expand Down Expand Up @@ -149,7 +150,7 @@ def determine_heat_flux(V_g, u, k, flux):

# Create output file
file_out = File("output/%s.pvd" % precice.get_participant_name())
file_out << u_n
file_out << (u_n, t)

print("output vtk for time = {}".format(float(t)))
n = 0
Expand All @@ -159,16 +160,16 @@ def determine_heat_flux(V_g, u, k, flux):

while precice.is_coupling_ongoing():

if precice.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint
if precice.requires_writing_checkpoint(): # write checkpoint
precice.store_checkpoint(u_n, t, n)

read_data = precice.read_data()
precice_dt = precice.get_max_time_step_size()
dt.assign(np.min([fenics_dt, precice_dt]))
read_data = precice.read_data(dt)

# Update the coupling expression with the new read data
precice.update_coupling_expression(coupling_expression, read_data)

dt.assign(np.min([fenics_dt, precice_dt]))

# Compute solution
solve(a == L, u_np1, bcs)

Expand All @@ -177,9 +178,9 @@ def determine_heat_flux(V_g, u, k, flux):
fluxes_y = fluxes.sub(1) # only exchange y component of flux.
precice.write_data(fluxes_y)

precice_dt = precice.advance(dt(0))
precice.advance(dt(0))

if precice.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint
if precice.requires_reading_checkpoint(): # roll back to checkpoint
u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
u_n.assign(u_cp)
t = t_cp
Expand All @@ -193,7 +194,7 @@ def determine_heat_flux(V_g, u, k, flux):
tol = 10e-5 # we need some tolerance, since otherwise output might be skipped.
if abs((t + tol) % dt_out) < 2 * tol: # output if t is a multiple of dt_out
print("output vtk for time = {}".format(float(t)))
file_out << u_n
file_out << (u_n, t)

# Update dirichlet BC
u_D.t = t + float(dt)
Expand Down
21 changes: 11 additions & 10 deletions partitioned-heat-conduction-complex/fenics/heat.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,14 +103,15 @@ def determine_gradient(V_g, u, flux):
# Initialize the adapter according to the specific participant
if problem is ProblemType.DIRICHLET:
precice = Adapter(adapter_config_filename="precice-adapter-config-D.json")
precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function)
precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function)
elif problem is ProblemType.NEUMANN:
precice = Adapter(adapter_config_filename="precice-adapter-config-N.json")
precice_dt = precice.initialize(coupling_boundary, read_function_space=V_g, write_object=u_D_function)
precice.initialize(coupling_boundary, read_function_space=V_g, write_object=u_D_function)

boundary_marker = False

dt = Constant(0)
precice_dt = precice.get_max_time_step_size()
dt.assign(np.min([fenics_dt, precice_dt]))

# Define variational problem
Expand Down Expand Up @@ -173,7 +174,7 @@ def determine_gradient(V_g, u, flux):
# 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
temperature_out << (u_n, t)
ref_out << u_ref
ranks << mesh_rank

Expand All @@ -191,10 +192,12 @@ def determine_gradient(V_g, u, flux):
while precice.is_coupling_ongoing():

# write checkpoint
if precice.is_action_required(precice.action_write_iteration_checkpoint()):
if precice.requires_writing_checkpoint():
precice.store_checkpoint(u_n, t, n)

read_data = precice.read_data()
precice_dt = precice.get_max_time_step_size()
dt.assign(np.min([fenics_dt, precice_dt]))
read_data = precice.read_data(dt)
if problem is ProblemType.DIRICHLET and (domain_part is DomainPart.CIRCULAR or domain_part is DomainPart.RECTANGLE):
# We have to data for an arbitrary point that is not on the circle, to obtain exact solution.
# See https://github.com/precice/fenics-adapter/issues/113 for details.
Expand All @@ -203,8 +206,6 @@ def determine_gradient(V_g, u, flux):
# Update the coupling expression with the new read data
precice.update_coupling_expression(coupling_expression, read_data)

dt.assign(np.min([fenics_dt, precice_dt]))

# Compute solution u^n+1, use bcs u_D^n+1, u^n and coupling bcs
solve(a == L, u_np1, bcs)

Expand All @@ -217,10 +218,10 @@ def determine_gradient(V_g, u, flux):
# Neumann problem reads flux and writes temperature on boundary to Dirichlet problem
precice.write_data(u_np1)

precice_dt = precice.advance(dt(0))
precice.advance(dt(0))

# roll back to checkpoint
if precice.is_action_required(precice.action_read_iteration_checkpoint()):
if precice.requires_reading_checkpoint():
u_cp, t_cp, n_cp = precice.retrieve_checkpoint()
u_n.assign(u_cp)
t = t_cp
Expand All @@ -237,7 +238,7 @@ def determine_gradient(V_g, u, flux):
print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error))
# output solution and reference solution at t_n+1
print('output u^%d and u_ref^%d' % (n, n))
temperature_out << u_n
temperature_out << (u_n, t)
ref_out << u_ref
error_out << error_pointwise

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
"config_file_name": "../precice-config.xml",
"interface": {
"coupling_mesh_name": "Dirichlet-Mesh",
"write_data_name": "Flux",
"write_data_name": "Heat-Flux",
"read_data_name": "Temperature"
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,6 @@
"interface": {
"coupling_mesh_name": "Neumann-Mesh",
"write_data_name": "Temperature",
"read_data_name": "Flux"
"read_data_name": "Heat-Flux"
}
}
14 changes: 7 additions & 7 deletions partitioned-heat-conduction-complex/precice-config.xml
Original file line number Diff line number Diff line change
Expand Up @@ -6,22 +6,22 @@
</log>

<data:scalar name="Temperature"/>
<data:vector name="Flux"/>
<data:vector name="Heat-Flux"/>

<mesh name="Dirichlet-Mesh" dimensions="2">
<use-data name="Temperature"/>
<use-data name="Flux"/>
<use-data name="Heat-Flux"/>
</mesh>

<mesh name="Neumann-Mesh" dimensions="2">
<use-data name="Temperature"/>
<use-data name="Flux"/>
<use-data name="Heat-Flux"/>
</mesh>

<participant name="Dirichlet">
<provide-mesh name="Dirichlet-Mesh" />
<receive-mesh name="Neumann-Mesh" from="Neumann" />
<write-data name="Flux" mesh="Dirichlet-Mesh" />
<write-data name="Heat-Flux" mesh="Dirichlet-Mesh" />
<read-data name="Temperature" mesh="Dirichlet-Mesh" />
<mapping:nearest-projection direction="read" from="Neumann-Mesh" to="Dirichlet-Mesh" constraint="consistent" />
</participant>
Expand All @@ -30,7 +30,7 @@
<provide-mesh name="Neumann-Mesh" />
<receive-mesh name="Dirichlet-Mesh" from="Dirichlet"/>
<write-data name="Temperature" mesh="Neumann-Mesh"/>
<read-data name="Flux" mesh="Neumann-Mesh"/>
<read-data name="Heat-Flux" mesh="Neumann-Mesh"/>
<mapping:nearest-projection direction="read" from="Dirichlet-Mesh" to="Neumann-Mesh" constraint="consistent" />
</participant>

Expand All @@ -41,9 +41,9 @@
<max-time value="1.0"/>
<time-window-size value="0.1"/>
<max-iterations value="100"/>
<exchange data="Flux" mesh="Dirichlet-Mesh" from="Dirichlet" to="Neumann" />
<exchange data="Heat-Flux" mesh="Dirichlet-Mesh" from="Dirichlet" to="Neumann" />
<exchange data="Temperature" mesh="Neumann-Mesh" from="Neumann" to="Dirichlet" initialize="true"/>
<relative-convergence-measure data="Flux" mesh="Dirichlet-Mesh" limit="1e-5"/>
<relative-convergence-measure data="Heat-Flux" mesh="Dirichlet-Mesh" limit="1e-5"/>
<relative-convergence-measure data="Temperature" mesh="Neumann-Mesh" limit="1e-5"/>
<acceleration:IQN-ILS>
<data name="Temperature" mesh="Neumann-Mesh"/>
Expand Down
Loading

0 comments on commit 63273dd

Please sign in to comment.