diff --git a/Notebooks/Examples-Meshing/Batmesh_number_one.h5 b/Notebooks/Examples-Meshing/Batmesh_number_one.h5 index c40eea8..ffa69db 100644 Binary files a/Notebooks/Examples-Meshing/Batmesh_number_one.h5 and b/Notebooks/Examples-Meshing/Batmesh_number_one.h5 differ diff --git a/Notebooks/Examples-Meshing/Convection-AdaptedMesh.py b/Notebooks/Examples-Meshing/Convection-AdaptedMesh.py index 7521ae6..c4302ca 100644 --- a/Notebooks/Examples-Meshing/Convection-AdaptedMesh.py +++ b/Notebooks/Examples-Meshing/Convection-AdaptedMesh.py @@ -39,12 +39,14 @@ class bd(Enum): batmesh.dm.view() # %% -Rayleigh_number = uw.function.expression(r"\textrm{Ra}", sympy.sympify(10)**6, "Rayleigh number") +Rayleigh_number = uw.function.expression(r"\textrm{Ra}", sympy.sympify(10000000), "Rayleigh number") # %% -v_soln = uw.discretisation.MeshVariable("U", batmesh, batmesh.dim, degree=2, continuous=True) -p_soln = uw.discretisation.MeshVariable("P", batmesh, 1, degree=1, continuous=True) -t_soln = uw.discretisation.MeshVariable("T", batmesh, 1, degree=2) +v_soln = uw.discretisation.MeshVariable("U", batmesh, batmesh.dim, degree=2, continuous=True) +v_soln1 = uw.discretisation.MeshVariable("U1", batmesh, batmesh.dim, degree=2, continuous=True) +p_soln = uw.discretisation.MeshVariable("P", batmesh, 1, degree=1, continuous=True) +t_soln = uw.discretisation.MeshVariable("T", batmesh, 1, degree=2, varsymbol=r"T_0") +t_soln1 = uw.discretisation.MeshVariable("T1", batmesh, 1, degree=2, varsymbol=r"T_1") cell_properties = uw.discretisation.MeshVariable("Bat", batmesh, 1, degree=0) @@ -68,18 +70,24 @@ class bd(Enum): stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 -stokes.tolerance = 1.0e-3 +stokes.tolerance = 1.0e-5 unit_r_vec = batmesh.CoordinateSystem.unit_e_0 +unit_y_vec = batmesh.CoordinateSystem.unit_j # free slip. # note with petsc we always need to provide a vector of correct cardinality. -stokes.add_natural_bc(10000 * unit_r_vec.dot(v_soln.sym) * unit_r_vec, "Upper") -stokes.bodyforce = Rayleigh_number * unit_r_vec * t_soln.sym[0] +stokes.add_essential_bc([0,0], "Upper") +# stokes.add_natural_bc(1e6 * Rayleigh_number * unit_r_vec.dot(v_soln.sym) * unit_r_vec, "Upper") -# Make the bat stagnant -# stokes.bodyforce -= 10000 * (1-cell_properties.sym[0]) * v_soln.sym +stokes.bodyforce = Rayleigh_number * unit_y_vec * t_soln.sym[0] + \ + Rayleigh_number * unit_y_vec * t_soln1.sym[0] + + +stokes.petsc_options.setValue("ksp_monitor", None) +stokes.petsc_options.setValue("snes_monitor", None) +stokes.petsc_options.setValue("snes_min_it", 1) stokes.view() @@ -92,31 +100,47 @@ class bd(Enum): adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel adv_diff.constitutive_model.Parameters.diffusivity = 1 -adv_diff.f = 10 * (1-cell_properties.sym[0]) - adv_diff.add_dirichlet_bc(0.0, "Upper") + +adv_diff1 = uw.systems.AdvDiffusionSLCN( + batmesh, + u_Field=t_soln1, + V_fn=v_soln1, +) + +adv_diff1.constitutive_model = uw.constitutive_models.DiffusionModel +adv_diff1.constitutive_model.Parameters.diffusivity = 1 +adv_diff1.add_dirichlet_bc(0.0, "Upper") + t_init = 1-cell_properties.sym[0] -with batmesh.access(t_soln): +with batmesh.access(t_soln, t_soln1): t_soln.data[:,0] = uw.function.evalf(t_init, t_soln.coords) + t_soln1.data[:,0] = uw.function.evalf(t_init, t_soln1.coords) -# %% -stokes.solve(zero_init_guess=True) # %% -v_soln.view() +stokes.solve(zero_init_guess=True, picard=0, verbose=False) +with batmesh.access(v_soln1): + v_soln1.data[...] = -v_soln.data[...] # %% + t_step = 0 # %% -for step in range(0, 50): +for step in range(0, 250): stokes.solve(zero_init_guess=False) - delta_t = 2.0 * stokes.estimate_dt() - adv_diff.solve(timestep=delta_t, zero_init_guess=False) + with batmesh.access(v_soln1): + v_soln1.data[...] = -v_soln.data[...] + + + delta_t = 5.0 * stokes.estimate_dt() + adv_diff.solve(timestep=delta_t, zero_init_guess=True) + adv_diff1.solve(timestep=delta_t, zero_init_guess=True) # stats then loop tstats = t_soln.stats() @@ -124,16 +148,17 @@ class bd(Enum): if uw.mpi.rank == 0: print("Timestep {}, dt {}".format(step, delta_t)) - batmesh.write_timestep( - "Batmesh", - meshUpdates=True, - meshVars=[p_soln, v_soln, t_soln], - outputPath="output", - index=t_step, - ) + if t_step%1==0: + batmesh.write_timestep( + "Batmesh", + meshUpdates=True, + meshVars=[p_soln, v_soln, t_soln, t_soln1], + outputPath="output", + index=t_step, + ) - with batmesh.access(t_soln): - t_soln.data[:,0] = np.maximum(t_soln.data[:,0], uw.function.evalf(t_init, t_soln.coords)) + # with batmesh.access(t_soln): + # t_soln.data[:,0] = np.maximum(t_soln.data[:,0], uw.function.evalf(t_init, t_soln.coords)) t_step += 1 @@ -146,13 +171,29 @@ class bd(Enum): pvmesh = vis.mesh_to_pv_mesh(batmesh) pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln.sym[0]) + pvmesh.point_data["R"] = vis.scalar_fn_to_pv_points(pvmesh, batmesh.CoordinateSystem.R[0]) + pvmesh.point_data["T1"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln1.sym[0]) pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym) pvmesh.point_data["L"] = vis.scalar_fn_to_pv_points(pvmesh, cell_properties.sym[0]) + pvmesh_t = vis.meshVariable_to_pv_mesh_object(t_soln) + pvmesh_t.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh_t, t_soln.sym[0]) + pvmesh_t.point_data["T1"] = vis.scalar_fn_to_pv_points(pvmesh_t, t_soln1.sym[0]) + + pvmesh_t.point_data["To"] = (pvmesh_t.point_data["T"] > 0.25).astype(float) + pvmesh_t.point_data["T1o"] = (pvmesh_t.point_data["T1"] > 0.25).astype(float) + pvmesh_t.point_data["TD"] = pvmesh_t.point_data["T"] - pvmesh_t.point_data["T1"] + + pl = pv.Plotter(window_size=(750, 750)) + pl.enable_depth_peeling() + doughnut = pvmesh.clip_scalar(scalars="R", + value=0.3, + invert=False) + pvstream = pvmesh.streamlines_from_source( - pvmesh.cell_centers(), + pv.PointSet(doughnut.cell_centers().points[::5]), vectors="V", integrator_type=45, surface_streamlines=True, @@ -160,55 +201,103 @@ class bd(Enum): max_time=0.25, ) + + # pl.add_mesh(pvmesh, + # style="wireframe", + # color="#FEFEF0", + # opacity=0.1) + + pl.add_mesh(pvstream, + cmap=["orange", "white", "green"], + show_scalar_bar=False, + opacity=0.1 + ) + + pl.add_mesh( - pvmesh, + pvmesh_t, cmap="Oranges", scalars="T", + opacity="sigmoid", edge_color="Grey", show_edges=False, use_transparency=False, - opacity=1.0, - show_scalar_bar=True, - clim=[0,1] + show_scalar_bar=False, ) pl.add_mesh( - pvmesh, - scalars="L", - cmap="Greys_r", - edge_color="Black", - opacity="L", + pvmesh_t, + copy_mesh=True, + cmap="Blues", + scalars="T1", + opacity="sigmoid", + edge_color="Grey", show_edges=False, - use_transparency=True, + use_transparency=False, show_scalar_bar=False, - + # clim=[0,2] ) - - # pl.add_mesh(pvstream, cmap="jet", show_scalar_bar=False) + + # pl.add_mesh( + # pvmesh, + # scalars="L", + # cmap="Greys_r", + # edge_color="Black", + # opacity="L", + # show_edges=False, + # use_transparency=True, + # show_scalar_bar=False, + # ) + # pl.add_arrows(pvmesh.points, pvmesh.point_data["V"], - # mag=float(20.0/Rayleigh_number.sym), + # mag=float(0.2/Rayleigh_number.sym), # show_scalar_bar=False) pl.show(jupyter_backend='html') +# %% +0/0 + +# %% + # %% ## Check the results saved to file expt_name = "Batmesh" output_path="output" -step = 20 +step = 50 v_soln.read_timestep(expt_name, "U", step, outputPath=output_path) p_soln.read_timestep(expt_name, "P", step, outputPath=output_path) t_soln.read_timestep(expt_name, "T", step, outputPath=output_path) +t_soln1.read_timestep(expt_name, "T1", step, outputPath=output_path) + +# %% +t_soln.stats() + +# %% +t_soln1.stats() # %% adv_diff # %% +adv_diff1 + +# %% +adv_diff1.DuDt.bdf(1) + +# %% +adv_diff.DuDt.bdf(1) + +# %% +batmesh.dm.view() + +# %% [markdown] +# # diff --git a/Notebooks/Examples-Meshing/MeshRefine-AdaptiveMetric-ImageDensity.py b/Notebooks/Examples-Meshing/MeshRefine-AdaptiveMetric-ImageDensity.py index 9561124..08db976 100644 --- a/Notebooks/Examples-Meshing/MeshRefine-AdaptiveMetric-ImageDensity.py +++ b/Notebooks/Examples-Meshing/MeshRefine-AdaptiveMetric-ImageDensity.py @@ -96,7 +96,7 @@ # %% with mesh0.access(H): - H.data[:,0] = 100.0 + 25000 * mvals + H.data[:,0] = 1000.0 + 50000 * mvals # + icoord, meshA = adaptivity.mesh_adapt_meshVar(mesh0, H, Metric, redistribute=True) @@ -129,97 +129,6 @@ # %% # ls -trl -# %% - -# %% - -# %% - -# %% -Rayleigh_number = uw.function.expression(r"\textrm{Ra}", sympy.sympify(10)**6, "Rayleigh number") - -# %% -v_soln = uw.discretisation.MeshVariable("U", batmesh, meshA.dim, degree=2, continuous=True) -p_soln = uw.discretisation.MeshVariable("P", batmesh, 1, degree=1, continuous=True) -t_soln = uw.discretisation.MeshVariable("T", batmesh, 1, degree=2) - - -# %% -# Create Stokes object - -stokes = uw.systems.Stokes( - batmesh, - velocityField=v_soln, - pressureField=p_soln, -) - -# Constant viscosity - -stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel -stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 -stokes.tolerance = 1.0e-3 - -unit_r_vec = batmesh.CoordinateSystem.unit_e_0 - -# free slip. -# note with petsc we always need to provide a vector of correct cardinality. - -stokes.add_natural_bc(10000 * unit_r_vec.dot(v_soln.sym) * unit_r_vec, "Upper") -stokes.bodyforce = -Rayleigh_number * unit_r_vec * t_soln.sym[0] - -stokes.view() - -# %% - -adv_diff = uw.systems.AdvDiffusionSLCN( - batmesh, - u_Field=t_soln, - V_fn=v_soln, -) - -adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel -adv_diff.constitutive_model.Parameters.diffusivity = 1 -adv_diff.f = cell_label.sym[0] - -adv_diff.add_dirichlet_bc(0.0, "Upper") - -with batmesh.access(t_soln): - t_soln.data[:,0] = uw.function.evalf(cell_label.sym[0], t_soln.coords) - - -# %% -stokes.solve(zero_init_guess=True) - -# %% -t_step = 0 - -# %% - -for step in range(0, 50): - stokes.solve(zero_init_guess=False) - delta_t = 5.0 * stokes.estimate_dt() - adv_diff.solve(timestep=delta_t, zero_init_guess=False) - - # stats then loop - tstats = t_soln.stats() - - if uw.mpi.rank == 0: - print("Timestep {}, dt {}".format(step, delta_t)) - - # if t_step % 5 == 0: - # plot_T_mesh(filename="{}_step_{}".format(expt_name, t_step)) - - batmesh.write_timestep( - "Batmesh", - meshUpdates=True, - meshVars=[p_soln, v_soln, t_soln], - outputPath="output", - index=t_step, - ) - - t_step += 1 - - # %% if uw.mpi.size == 1: @@ -227,24 +136,14 @@ import underworld3.visualisation as vis pvmesh = vis.mesh_to_pv_mesh(batmesh) - pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln.sym[0]) - pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym) + pvmesh.point_data["L"] = vis.scalar_fn_to_pv_points(pvmesh, cell_label.sym[0]) pl = pv.Plotter(window_size=(1500, 1500)) - pvstream = pvmesh.streamlines_from_source( - pvmesh.cell_centers(), - vectors="V", - integrator_type=45, - surface_streamlines=True, - max_steps=1000, - max_time=0.25, - ) - pl.add_mesh( pvmesh, cmap="Oranges_r", - scalars="T", + scalars="L", edge_color="Black", show_edges=True, use_transparency=False, @@ -264,9 +163,7 @@ ) - pl.add_mesh(pvstream, cmap="Oranges") - pl.add_arrows(pvmesh.points, pvmesh.point_data["V"], mag=float(10.0/Rayleigh_number.sym)) pl.show(jupyter_backend='html')