From f693c9fa73bbb713c0a0f1644adec9faa520b7c9 Mon Sep 17 00:00:00 2001 From: Henrik Kjedsberg Date: Thu, 4 Jan 2024 16:58:03 +0100 Subject: [PATCH 1/3] Add cylinder mesh and use mesh_path for Cylinder problem (rigid) --- src/oasismove/problems/Cylinder.py | 31 +++++-------------- src/oasismove/problems/NSCoupled/Cylinder.py | 17 ++++++++-- src/oasismove/problems/NSfracStep/Cylinder.py | 18 ++++++++--- 3 files changed, 35 insertions(+), 31 deletions(-) diff --git a/src/oasismove/problems/Cylinder.py b/src/oasismove/problems/Cylinder.py index 990cf7e..99175ce 100644 --- a/src/oasismove/problems/Cylinder.py +++ b/src/oasismove/problems/Cylinder.py @@ -3,28 +3,7 @@ __copyright__ = "Copyright (C) 2014 " + __author__ __license__ = "GNU Lesser GPL version 3 or any later version" -import os -import platform - -from dolfin import Mesh, AutoSubDomain, near - -if not os.path.isfile("cylinder.xml"): - if platform.system() == "Linux": - os.system("wget -O cylinder.xml https://www.dropbox.com/s/d78g4cyjxl3ylay/cylinder.xml?dl=0") - elif platform.system() == "Darwin": - os.system("curl -L https://www.dropbox.com/s/d78g4cyjxl3ylay/cylinder.xml?dl=0 -o cylinder.xml") - else: - raise ImportError("Could not determine platform") - - # try: - # os.system("gmsh mesh/cylinder.geo -2 -o mesh/cylinder.msh") - # os.system("dolfin-convert mesh/cylinder.msh mesh/cylinder.xml") - # os.system("rm mesh/cylinder.msh") - # except RuntimeError: - # os.system("wget -O cylinder.xml https://www.dropbox.com/s/d78g4cyjxl3ylay/cylinder.xml?dl=0") - # raise "Gmsh is required to run this demo" - -mesh = Mesh("cylinder.xml") +from dolfin import AutoSubDomain, near H = 0.41 L = 2.2 @@ -47,7 +26,11 @@ # Overload post_import_problem to choose between the two cases -def post_import_problem(NS_parameters, commandline_kwargs, **NS_namespace): +def post_import_problem(NS_parameters, mesh, commandline_kwargs, **NS_namespace): + # If the mesh is a callable function, then create the mesh here. + if callable(mesh): + mesh = mesh(**NS_parameters) + """ Choose case - case could be defined through command line.""" NS_parameters.update(commandline_kwargs) case = NS_parameters['case'] if 'case' in NS_parameters else 1 @@ -55,6 +38,6 @@ def post_import_problem(NS_parameters, commandline_kwargs, **NS_namespace): Re = cases[case]["Re"] Umean = 2. / 3. * Um nu = Umean * D / Re - NS_parameters.update(nu=nu, Re=Re, Um=Um, Umean=Umean) + NS_parameters.update(nu=nu, Re=Re, Um=Um, Umean=Umean, mesh=mesh) return NS_parameters diff --git a/src/oasismove/problems/NSCoupled/Cylinder.py b/src/oasismove/problems/NSCoupled/Cylinder.py index d2094c4..c02639c 100644 --- a/src/oasismove/problems/NSCoupled/Cylinder.py +++ b/src/oasismove/problems/NSCoupled/Cylinder.py @@ -5,8 +5,8 @@ __copyright__ = "Copyright (C) 2014 " + __author__ __license__ = "GNU Lesser GPL version 3 or any later version" -from ..Cylinder import * -from ..NSCoupled import * +from oasismove.problems.NSCoupled import * +from oasismove.problems.Cylinder import * # Override some problem specific parameters @@ -15,7 +15,9 @@ def problem_parameters(NS_parameters, scalar_components, **NS_namespace): omega=1.0, max_iter=100, plot_interval=10, - velocity_degree=2) + velocity_degree=2, + mesh_path="src/oasismove/mesh/cylinder.xml", + ) scalar_components += ["c", "d"] @@ -24,6 +26,15 @@ def scalar_source(c_, d_, **NS_namespace): return {"c": -Constant(0.1) * c_ * c_, "d": -Constant(0.25) * c_ * d_ * d_} +def mesh(mesh_path, dt, **NS_namespace): + # Import mesh + print(mesh_path) + mesh = Mesh(mesh_path) + + print_mesh_information(mesh, dt, dim=2) + return mesh + + def create_bcs(VQ, Um, CG, V, element, **NS_namespace): inlet = Expression(("4.*{0}*x[1]*({1}-x[1])/pow({1}, 2)".format(Um, H), "0"), element=V) ux = Expression(("0.00*x[1]", "-0.00*(x[0]-{})".format(center)), element=V) diff --git a/src/oasismove/problems/NSfracStep/Cylinder.py b/src/oasismove/problems/NSfracStep/Cylinder.py index 25f87fa..8a2e912 100644 --- a/src/oasismove/problems/NSfracStep/Cylinder.py +++ b/src/oasismove/problems/NSfracStep/Cylinder.py @@ -9,9 +9,8 @@ from os import getcwd import matplotlib.pyplot as plt - -from ..Cylinder import * -from ..NSfracStep import * +from oasismove.problems.NSfracStep import * +from oasismove.problems.Cylinder import * def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, @@ -35,7 +34,9 @@ def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, plot_interval=10, velocity_degree=2, print_intermediate_info=100, - use_krylov_solvers=True) + use_krylov_solvers=True, + mesh_path="src/oasismove/mesh/cylinder.xml", + ) NS_parameters['krylov_solvers'].update(dict(monitor_convergence=True)) NS_parameters['velocity_krylov_solver'].update(dict(preconditioner_type='jacobi', solver_type='bicgstab')) @@ -44,6 +45,15 @@ def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, Schmidt["alfa"] = 0.1 +def mesh(mesh_path, dt, **NS_namespace): + # Import mesh + print(mesh_path) + mesh = Mesh(mesh_path) + + print_mesh_information(mesh, dt, dim=2) + return mesh + + def create_bcs(V, Q, Um, H, **NS_namespace): inlet = Expression( "4.*{0}*x[1]*({1}-x[1])/pow({1}, 2)".format(Um, H), degree=2) From d00367a4c64e4cb41cea71f35270eb26db71a5f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Mon, 8 Jan 2024 11:19:36 +0100 Subject: [PATCH 2/3] Fix path input for cylinder example --- src/oasismove/NSfracStep.py | 3 ++- src/oasismove/common/io.py | 1 + src/oasismove/problems/NSCoupled/Cylinder.py | 13 ++++++++++--- src/oasismove/problems/NSfracStep/Cylinder.py | 11 ++++++++--- 4 files changed, 21 insertions(+), 7 deletions(-) diff --git a/src/oasismove/NSfracStep.py b/src/oasismove/NSfracStep.py index 7dec941..3d5f1a8 100755 --- a/src/oasismove/NSfracStep.py +++ b/src/oasismove/NSfracStep.py @@ -48,7 +48,7 @@ problemspec = importlib.util.find_spec(problemname) if problemspec is None: raise RuntimeError(problemname + ' not found') - +mesh_path = commandline_kwargs.get('mesh_path', None) # Import the problem module print('Importing problem module ' + problemname + ':\n' + problemspec.origin) problemmod = importlib.util.module_from_spec(problemspec) @@ -77,6 +77,7 @@ tstep = params["tstep"] previous_velocity_degree = params["velocity_degree"] + # Import chosen functionality from solvers solver = importlib.import_module('.'.join(('oasismove.solvers.NSfracStep', solver))) vars().update({name: solver.__dict__[name] for name in solver.__all__}) diff --git a/src/oasismove/common/io.py b/src/oasismove/common/io.py index a00e317..662ab6f 100644 --- a/src/oasismove/common/io.py +++ b/src/oasismove/common/io.py @@ -114,6 +114,7 @@ def save_tstep_solution_xdmf(tstep, q_, u_, newfolder, tstepfiles, output_timese if MPI.rank(MPI.comm_world) == 0: if not path.exists(path.join(timefolder, "params.dat")): f = open(path.join(timefolder, 'params.dat'), 'wb') + NS_parameters.pop("mesh") pickle.dump(NS_parameters, f) diff --git a/src/oasismove/problems/NSCoupled/Cylinder.py b/src/oasismove/problems/NSCoupled/Cylinder.py index c02639c..6c8938a 100644 --- a/src/oasismove/problems/NSCoupled/Cylinder.py +++ b/src/oasismove/problems/NSCoupled/Cylinder.py @@ -7,16 +7,18 @@ from oasismove.problems.NSCoupled import * from oasismove.problems.Cylinder import * - +from pathlib import Path # Override some problem specific parameters + + def problem_parameters(NS_parameters, scalar_components, **NS_namespace): NS_parameters.update( omega=1.0, max_iter=100, plot_interval=10, velocity_degree=2, - mesh_path="src/oasismove/mesh/cylinder.xml", + mesh_path="src/oasismove/mesh/cylinder.xdmf", ) scalar_components += ["c", "d"] @@ -29,7 +31,12 @@ def scalar_source(c_, d_, **NS_namespace): def mesh(mesh_path, dt, **NS_namespace): # Import mesh print(mesh_path) - mesh = Mesh(mesh_path) + if Path(mesh_path).suffix == ".xml": + mesh = Mesh(mesh_path) + elif Path(mesh_path).suffix == ".xdmf": + mesh = Mesh() + with XDMFFile(mesh_path) as infile: + infile.read(mesh) print_mesh_information(mesh, dt, dim=2) return mesh diff --git a/src/oasismove/problems/NSfracStep/Cylinder.py b/src/oasismove/problems/NSfracStep/Cylinder.py index 8a2e912..79bde06 100644 --- a/src/oasismove/problems/NSfracStep/Cylinder.py +++ b/src/oasismove/problems/NSfracStep/Cylinder.py @@ -11,6 +11,7 @@ import matplotlib.pyplot as plt from oasismove.problems.NSfracStep import * from oasismove.problems.Cylinder import * +from pathlib import Path def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, @@ -23,7 +24,6 @@ def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, NS_parameters.update(pickle.load(f)) NS_parameters['restart_folder'] = restart_folder globals().update(NS_parameters) - else: # Override some problem specific parameters NS_parameters.update( @@ -35,7 +35,7 @@ def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, velocity_degree=2, print_intermediate_info=100, use_krylov_solvers=True, - mesh_path="src/oasismove/mesh/cylinder.xml", + mesh_path=commandline_kwargs["mesh_path"], ) NS_parameters['krylov_solvers'].update(dict(monitor_convergence=True)) NS_parameters['velocity_krylov_solver'].update(dict(preconditioner_type='jacobi', @@ -48,7 +48,12 @@ def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, def mesh(mesh_path, dt, **NS_namespace): # Import mesh print(mesh_path) - mesh = Mesh(mesh_path) + if Path(mesh_path).suffix == ".xml": + mesh = Mesh(mesh_path) + elif Path(mesh_path).suffix == ".xdmf": + mesh = Mesh() + with XDMFFile(mesh_path) as infile: + infile.read(mesh) print_mesh_information(mesh, dt, dim=2) return mesh From ba44e419d857a9cfd3f054c8c62f83588d1886b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Mon, 8 Jan 2024 11:21:18 +0100 Subject: [PATCH 3/3] Remove mesh path from file --- src/oasismove/NSfracStep.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oasismove/NSfracStep.py b/src/oasismove/NSfracStep.py index 3d5f1a8..3119942 100755 --- a/src/oasismove/NSfracStep.py +++ b/src/oasismove/NSfracStep.py @@ -48,7 +48,7 @@ problemspec = importlib.util.find_spec(problemname) if problemspec is None: raise RuntimeError(problemname + ' not found') -mesh_path = commandline_kwargs.get('mesh_path', None) + # Import the problem module print('Importing problem module ' + problemname + ':\n' + problemspec.origin) problemmod = importlib.util.module_from_spec(problemspec)