diff --git a/partitioned-heat-conduction/doConvergenceStudy.py b/partitioned-heat-conduction/doConvergenceStudy.py new file mode 100644 index 000000000..e8e8e93ee --- /dev/null +++ b/partitioned-heat-conduction/doConvergenceStudy.py @@ -0,0 +1,81 @@ +from jinja2 import Environment, select_autoescape, FileSystemLoader +import pandas as pd +from pathlib import Path +import subprocess +import datetime +import os + + +def render(dt): + base_path = Path(__file__).parent.absolute() + + env = Environment( + loader=FileSystemLoader(base_path), + autoescape=select_autoescape(['xml']) + ) + + precice_config_template = env.get_template('precice-config-template.xml') + + precice_config_name = base_path / "precice-config.xml" + + with open(precice_config_name, "w") as file: + file.write(precice_config_template.render(time_window_size=dt)) + + +def do_run(dt, error_tol=10e-3): + fenics = Path(__file__).parent.absolute() / "fenics" + render(dt=dt) + print(f"Start run with dt={dt} at {datetime.datetime.now()} ...") + participants = [ + { + "name": "Dirichlet", + "cmd":"-d", + }, + { + "name": "Neumann", + "cmd":"-n", + }, + ] + + for participant in participants: + with open(fenics / f"stdout-{participant['name']}.log", "w") as outfile: + p = subprocess.Popen(["python3", fenics / "heat.py", participant["cmd"], f"-e {error_tol}"], cwd=fenics, stdout=outfile) + participant["proc"] = p + + for participant in participants: + participant["proc"].wait() + + for participant in participants: + if participant["proc"].returncode != 0: + raise Exception(f'Experiment with dt={dt} failed!') + + print("Done.") + print("Postprocessing...") + summary = {"dt":dt} + for participant in participants: + df = pd.read_csv(fenics / f"errors-{participant['name']}.csv", comment="#") + summary[f"error {participant['name']}"] = df["errors"].abs().max() + print("Done.") + + return summary + + +if __name__ == "__main__": + min_dt = 0.01 + dts = [min_dt * 0.5**i for i in range(10)] + + df = pd.DataFrame(columns=["dt", "error Dirichlet", "error Neumann"]) + + summary_file = "convergence_study.csv" + for dt in dts: + summary = do_run(dt, error_tol=10e10) + df = pd.concat([df, pd.DataFrame(summary, index=[0])], ignore_index=True) + + print(f"Write output to {summary_file}") + df.to_csv(summary_file) + + term_size = os.get_terminal_size() + print('-' * term_size.columns) + print("Preliminary results:") + print(df) + print('-' * term_size.columns) \ No newline at end of file diff --git a/partitioned-heat-conduction/fenics/heat.py b/partitioned-heat-conduction/fenics/heat.py index d3a35a2a2..45c566ae8 100644 --- a/partitioned-heat-conduction/fenics/heat.py +++ b/partitioned-heat-conduction/fenics/heat.py @@ -36,6 +36,7 @@ from problem_setup import get_geometry, get_manufactured_solution import dolfin from dolfin import FacetNormal, dot +import pandas as pd def determine_gradient(V_g, u, flux): @@ -109,6 +110,7 @@ def determine_gradient(V_g, u, flux): precice_dt = precice.initialize(coupling_boundary, read_function_space=W, write_object=u_D_function) dt = Constant(0) +window_dt = precice_dt # store for later fenics_dt = precice_dt # use window size provided by preCICE as time step size dt.assign(np.min([fenics_dt, precice_dt])) @@ -166,6 +168,8 @@ def determine_gradient(V_g, u, flux): error_total, error_pointwise = compute_errors(u_n, u_ref, V) error_out << error_pointwise +errors = [] +times = [] # set t_1 = t_0 + dt, this gives u_D^1 # call dt(0) to evaluate FEniCS Constant. Todo: is there a better way? @@ -225,6 +229,8 @@ def determine_gradient(V_g, u, flux): temperature_out << u_n ref_out << u_ref error_out << error_pointwise + errors.append(error) + times.append(t) # Update Dirichlet BC u_D.t = t + float(dt) @@ -232,3 +238,21 @@ def determine_gradient(V_g, u, flux): # Hold plot precice.finalize() + +df = pd.DataFrame() +df["times"] = times +df["errors"] = errors +df = df.set_index('times') +metadata = f'''# time_window_size: {window_dt} +# time_step_size: {fenics_dt} +''' + +filename = f"errors-{problem.value}.csv" + +import os + +os.remove(filename) + +with open(filename, 'a') as f: + f.write(f"{metadata}") + df.to_csv(f) \ No newline at end of file diff --git a/partitioned-heat-conduction/fenics/my_enums.py b/partitioned-heat-conduction/fenics/my_enums.py index 9b203aeb2..42120fc78 100644 --- a/partitioned-heat-conduction/fenics/my_enums.py +++ b/partitioned-heat-conduction/fenics/my_enums.py @@ -5,8 +5,8 @@ class ProblemType(Enum): """ Enum defines problem type. Details see above. """ - DIRICHLET = 1 # Dirichlet problem - NEUMANN = 2 # Neumann problem + DIRICHLET = "Dirichlet" # Dirichlet problem + NEUMANN = "Neumann" # Neumann problem class DomainPart(Enum): diff --git a/partitioned-heat-conduction/precice-config.xml b/partitioned-heat-conduction/precice-config-template.xml similarity index 86% rename from partitioned-heat-conduction/precice-config.xml rename to partitioned-heat-conduction/precice-config-template.xml index 280fcb2f2..a27fb0a2c 100644 --- a/partitioned-heat-conduction/precice-config.xml +++ b/partitioned-heat-conduction/precice-config-template.xml @@ -52,7 +52,7 @@ - + - - + + + - + +