diff --git a/src/oasismove/NSfracStep.py b/src/oasismove/NSfracStep.py index 7dec941..3119942 100755 --- a/src/oasismove/NSfracStep.py +++ b/src/oasismove/NSfracStep.py @@ -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/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..6c8938a 100644 --- a/src/oasismove/problems/NSCoupled/Cylinder.py +++ b/src/oasismove/problems/NSCoupled/Cylinder.py @@ -5,17 +5,21 @@ __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 * +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) + velocity_degree=2, + mesh_path="src/oasismove/mesh/cylinder.xdmf", + ) scalar_components += ["c", "d"] @@ -24,6 +28,20 @@ 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) + 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 + + 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..79bde06 100644 --- a/src/oasismove/problems/NSfracStep/Cylinder.py +++ b/src/oasismove/problems/NSfracStep/Cylinder.py @@ -9,9 +9,9 @@ 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 * +from pathlib import Path def problem_parameters(commandline_kwargs, NS_parameters, scalar_components, @@ -24,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 +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=commandline_kwargs["mesh_path"], + ) 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,20 @@ 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) + 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 + + 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)