-
-
Notifications
You must be signed in to change notification settings - Fork 109
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Restructure partitioned-heat-conduction-direct (#497)
Co-authored-by: Jun Chen <jun.chen@ipvs.uni-stuttgart.de>
- Loading branch information
Showing
11 changed files
with
314 additions
and
190 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
File renamed without changes.
151 changes: 151 additions & 0 deletions
151
partitioned-heat-conduction-direct/dirichlet-nutils/heat.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,151 @@ | ||
#! /usr/bin/env python3 | ||
|
||
from nutils import cli, mesh, function, solver, export | ||
import functools | ||
import treelog | ||
import numpy as np | ||
import precice | ||
|
||
|
||
def main(n=10, degree=1, timestep=.1, alpha=3., beta=1.2): | ||
|
||
x_grid = np.linspace(0, 1, n) | ||
y_grid = np.linspace(0, 1, n) | ||
|
||
# define the Nutils mesh | ||
domain, geom = mesh.rectilinear([x_grid, y_grid]) | ||
coupling_boundary = domain.boundary['right'] | ||
read_sample = coupling_boundary.sample('gauss', degree=degree * 2) | ||
|
||
# Nutils namespace | ||
ns = function.Namespace() | ||
ns.x = geom | ||
ns.basis = domain.basis('std', degree=degree) | ||
ns.alpha = alpha # parameter of problem | ||
ns.beta = beta # parameter of problem | ||
ns.u = 'basis_n ?lhs_n' # solution | ||
ns.dudt = 'basis_n (?lhs_n - ?lhs0_n) / ?dt' # time derivative | ||
ns.flux = 'basis_n ?fluxdofs_n' # heat flux | ||
ns.f = 'beta - 2 - 2 alpha' # rhs | ||
ns.uexact = '1 + x_0 x_0 + alpha x_1 x_1 + beta ?t' # analytical solution | ||
ns.readbasis = read_sample.basis() | ||
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) | ||
|
||
# set boundary conditions at non-coupling boundaries | ||
# top and bottom boundary are non-coupling for both sides | ||
sqr = domain.boundary['top,bottom,left'].integral( | ||
'(u - uexact)^2 d:x' @ ns, degree=degree * 2) | ||
|
||
sqr += read_sample.integral('(u - readfunc)^2 d:x' @ ns) | ||
|
||
# preCICE setup | ||
participant = precice.Participant( | ||
"Dirichlet", "../precice-config.xml", 0, 1) | ||
|
||
mesh_name_read = "Dirichlet-Mesh" | ||
mesh_name_write = "Neumann-Mesh" | ||
|
||
vertex_ids_read = participant.set_mesh_vertices( | ||
mesh_name_read, read_sample.eval(ns.x)) | ||
participant.set_mesh_access_region(mesh_name_write, [.9, 1.1, -.1, 1.1]) | ||
|
||
participant.initialize() | ||
precice_dt = participant.get_max_time_step_size() | ||
solver_dt = timestep | ||
dt = min(precice_dt, solver_dt) | ||
|
||
vertex_ids_write, coords = participant.get_mesh_vertex_ids_and_coordinates( | ||
mesh_name_write) | ||
write_sample = domain.locate(ns.x, coords, eps=1e-10, tol=1e-10) | ||
|
||
precice_write = functools.partial( | ||
participant.write_data, mesh_name_write, "Heat-Flux", vertex_ids_write) | ||
precice_read = functools.partial( | ||
participant.read_data, mesh_name_read, "Temperature", vertex_ids_read) | ||
|
||
# helper functions to project heat flux to coupling boundary | ||
|
||
# To communicate the flux to the Neumann side we should not simply | ||
# evaluate u_,i n_i as this is an unbounded term leading to suboptimal | ||
# convergence. Instead we project ∀ v: ∫_Γ v flux = ∫_Γ v u_,i n_i and | ||
# evaluate flux. While the right-hand-side contains the same unbounded | ||
# term, we can use the strong identity du/dt - u_,ii = f to rewrite it | ||
# to ∫_Ω [v (du/dt - f) + v_,i u_,i] - ∫_∂Ω\Γ v u_,k n_k, in which we | ||
# recognize the residual and an integral over the exterior boundary. | ||
# 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. | ||
right_sqr = domain.boundary['right'].integral( | ||
'flux^2 d:x' @ ns, degree=degree * 2) | ||
right_cons = solver.optimize('fluxdofs', right_sqr, droptol=1e-10) | ||
# right_cons is NaN in dofs that are NOT supported on the right boundary | ||
flux_sqr = domain.boundary['right'].boundary['top,bottom'].integral( | ||
'(flux - uexact_,0)^2 d:x' @ ns, degree=degree * 2) | ||
flux_cons = solver.optimize('fluxdofs', flux_sqr, droptol=1e-10, | ||
constrain=np.choose(np.isnan(right_cons), [np.nan, 0.])) | ||
# flux_cons is NaN in dofs that are supported on ONLY the right boundary | ||
flux_res = read_sample.integral('basis_n flux d:x' @ ns) - res | ||
|
||
t = 0. | ||
istep = 0 | ||
|
||
# initial condition | ||
sqr0 = domain.integral('(u - uexact)^2' @ ns, degree=degree * 2) | ||
lhs = solver.optimize('lhs', sqr0, arguments=dict(t=t)) | ||
bezier = domain.sample('bezier', degree * 2) | ||
|
||
while participant.is_coupling_ongoing(): | ||
|
||
# save checkpoint | ||
if participant.requires_writing_checkpoint(): | ||
checkpoint = lhs, t, istep | ||
|
||
# prepare next timestep | ||
precice_dt = participant.get_max_time_step_size() | ||
dt = min(precice_dt, solver_dt) | ||
lhs0 = lhs | ||
istep += 1 | ||
t += dt | ||
|
||
# read data from participant | ||
read_data = precice_read(dt) | ||
|
||
# update (time-dependent) boundary condition | ||
cons = solver.optimize('lhs', sqr, droptol=1e-15, | ||
arguments=dict(t=t, readdata=read_data)) | ||
|
||
# solve nutils timestep | ||
lhs = solver.solve_linear('lhs', res, constrain=cons, arguments=dict( | ||
lhs0=lhs0, dt=dt, t=t, readdata=read_data)) | ||
|
||
# write data to participant | ||
fluxdofs = solver.solve_linear( | ||
'fluxdofs', flux_res, arguments=dict( | ||
lhs0=lhs0, lhs=lhs, dt=dt, t=t), constrain=flux_cons) | ||
write_data = write_sample.eval('flux' @ ns, fluxdofs=fluxdofs) | ||
|
||
precice_write(write_data) | ||
|
||
# do the coupling | ||
participant.advance(dt) | ||
|
||
# read checkpoint if required | ||
if participant.requires_reading_checkpoint(): | ||
lhs, t, istep = checkpoint | ||
else: | ||
# generate output | ||
x, u, uexact = bezier.eval( | ||
['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t) | ||
with treelog.add(treelog.DataLog()): | ||
export.vtk("Dirichlet" + "-" + str(istep), bezier.tri, | ||
x, Temperature=u, reference=uexact) | ||
|
||
participant.finalize() | ||
|
||
|
||
if __name__ == '__main__': | ||
cli.run(main) |
File renamed without changes.
14 changes: 14 additions & 0 deletions
14
partitioned-heat-conduction-direct/dirichlet-nutils/run.sh
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,14 @@ | ||
#!/bin/bash | ||
set -e -u | ||
|
||
. ../../tools/log.sh | ||
exec > >(tee --append "$LOGFILE") 2>&1 | ||
|
||
python3 -m venv .venv | ||
. .venv/bin/activate | ||
pip install -r requirements.txt | ||
|
||
rm -rf Dirichlet-*.vtk | ||
NUTILS_RICHOUTPUT=no python3 heat.py | ||
|
||
close_log |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,6 @@ | ||
#!/bin/sh | ||
set -e -u | ||
|
||
. ../../tools/cleaning-tools.sh | ||
|
||
clean_nutils . |
123 changes: 123 additions & 0 deletions
123
partitioned-heat-conduction-direct/neumann-nutils/heat.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,123 @@ | ||
#! /usr/bin/env python3 | ||
|
||
from nutils import cli, mesh, function, solver, export | ||
import functools | ||
import treelog | ||
import numpy as np | ||
import precice | ||
|
||
|
||
def main(n=10, degree=1, timestep=.1, alpha=3., beta=1.2): | ||
|
||
x_grid = np.linspace(1, 2, n) | ||
y_grid = np.linspace(0, 1, n) | ||
|
||
# define the Nutils mesh | ||
domain, geom = mesh.rectilinear([x_grid, y_grid]) | ||
coupling_boundary = domain.boundary['left'] | ||
read_sample = coupling_boundary.sample('gauss', degree=degree * 2) | ||
|
||
# Nutils namespace | ||
ns = function.Namespace() | ||
ns.x = geom | ||
ns.basis = domain.basis('std', degree=degree) | ||
ns.alpha = alpha # parameter of problem | ||
ns.beta = beta # parameter of problem | ||
ns.u = 'basis_n ?lhs_n' # solution | ||
ns.dudt = 'basis_n (?lhs_n - ?lhs0_n) / ?dt' # time derivative | ||
ns.flux = 'basis_n ?fluxdofs_n' # heat flux | ||
ns.f = 'beta - 2 - 2 alpha' # rhs | ||
ns.uexact = '1 + x_0 x_0 + alpha x_1 x_1 + beta ?t' # analytical solution | ||
ns.readbasis = read_sample.basis() | ||
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) | ||
|
||
# set boundary conditions at non-coupling boundaries | ||
# top and bottom boundary are non-coupling for both sides | ||
sqr = domain.boundary['top,bottom,right'].integral( | ||
'(u - uexact)^2 d:x' @ ns, degree=degree * 2) | ||
|
||
res += read_sample.integral('basis_n readfunc d:x' @ ns) | ||
|
||
# preCICE setup | ||
participant = precice.Participant("Neumann", "../precice-config.xml", 0, 1) | ||
|
||
mesh_name_read = "Neumann-Mesh" | ||
mesh_name_write = "Dirichlet-Mesh" | ||
|
||
vertex_ids_read = participant.set_mesh_vertices( | ||
mesh_name_read, read_sample.eval(ns.x)) | ||
participant.set_mesh_access_region(mesh_name_write, [.9, 1.1, -.1, 1.1]) | ||
|
||
participant.initialize() | ||
precice_dt = participant.get_max_time_step_size() | ||
solver_dt = timestep | ||
dt = min(precice_dt, solver_dt) | ||
|
||
vertex_ids_write, coords = participant.get_mesh_vertex_ids_and_coordinates( | ||
mesh_name_write) | ||
write_sample = domain.locate(ns.x, coords, eps=1e-10, tol=1e-10) | ||
|
||
precice_write = functools.partial( | ||
participant.write_data, mesh_name_write, "Temperature", vertex_ids_write) | ||
precice_read = functools.partial( | ||
participant.read_data, mesh_name_read, "Heat-Flux", vertex_ids_read) | ||
|
||
t = 0. | ||
istep = 0 | ||
|
||
# initial condition | ||
sqr0 = domain.integral('(u - uexact)^2' @ ns, degree=degree * 2) | ||
lhs = solver.optimize('lhs', sqr0, arguments=dict(t=t)) | ||
bezier = domain.sample('bezier', degree * 2) | ||
|
||
while participant.is_coupling_ongoing(): | ||
|
||
# save checkpoint | ||
if participant.requires_writing_checkpoint(): | ||
checkpoint = lhs, t, istep | ||
|
||
# prepare next timestep | ||
precice_dt = participant.get_max_time_step_size() | ||
dt = min(precice_dt, solver_dt) | ||
lhs0 = lhs | ||
istep += 1 | ||
t += dt | ||
|
||
# read data from participant | ||
read_data = precice_read(dt) | ||
|
||
# update (time-dependent) boundary condition | ||
cons = solver.optimize('lhs', sqr, droptol=1e-15, | ||
arguments=dict(t=t, readdata=read_data)) | ||
|
||
# solve nutils timestep | ||
lhs = solver.solve_linear('lhs', res, constrain=cons, arguments=dict( | ||
lhs0=lhs0, dt=dt, t=t, readdata=read_data)) | ||
|
||
# write data to participant | ||
write_data = write_sample.eval('u' @ ns, lhs=lhs) | ||
precice_write(write_data) | ||
|
||
# do the coupling | ||
participant.advance(dt) | ||
|
||
# read checkpoint if required | ||
if participant.requires_reading_checkpoint(): | ||
lhs, t, istep = checkpoint | ||
else: | ||
# generate output | ||
x, u, uexact = bezier.eval( | ||
['x_i', 'u', 'uexact'] @ ns, lhs=lhs, t=t) | ||
with treelog.add(treelog.DataLog()): | ||
export.vtk("Neumann" + "-" + str(istep), bezier.tri, | ||
x, Temperature=u, reference=uexact) | ||
|
||
participant.finalize() | ||
|
||
|
||
if __name__ == '__main__': | ||
cli.run(main) |
2 changes: 2 additions & 0 deletions
2
partitioned-heat-conduction-direct/neumann-nutils/requirements.txt
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,2 @@ | ||
nutils==7 | ||
pyprecice==3 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,14 @@ | ||
#!/bin/bash | ||
set -e -u | ||
|
||
. ../../tools/log.sh | ||
exec > >(tee --append "$LOGFILE") 2>&1 | ||
|
||
python3 -m venv .venv | ||
. .venv/bin/activate | ||
pip install -r requirements.txt | ||
|
||
rm -rf Neumann-*.vtk | ||
NUTILS_RICHOUTPUT=no python3 heat.py | ||
|
||
close_log |
Oops, something went wrong.