omega_h
omega_h copied to clipboard
warning: gradation limiting is up to step 2200
This case is similar to the case described in question #244 but in 3d.
As one can see from the logfile, the first outer iteration (marked with "step: 1") runs through smoothly. In contrast, the second step takes very long (generated metrics in 0 iterations and 587.519 seconds
) and reports a lot of warning: gradation limiting is up to step ...
. Interestingly, steps 3 to 10 run through smoothly again.
@ibaned Do you have any pointers how to improve the code or is this the expected behaviour in 3d?
As only the second adaption step takes that long, it is not a major problem. For production runs I would take the following approach:
- Start with a very coarse and uniform mesh (mesh1)
- Interpolate the initial conditions on mesh1
- Adapt mesh1 and get mesh2
- Interpolate the initial conditions on mesh2
- Repeat steps 3 and 4 a few times
- Write final adapted mesh to disk and reuse this mesh in subsequent production runs.
@ibaned Do you think, that this is a sound approach?
log: demo_omegah_dolfin_3d.log
import PyOmega_h as omega_h
import dolfin as dolfin
PHASEFIELD = "phi"
# initial conditions
f_code = """
#include <iostream>
#include <cmath>
#include <pybind11/pybind11.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
#include <dolfin/function/Constant.h>
#include <dolfin/common/Array.h>
class F : public dolfin::Expression {
public:
// Create expression
F() : dolfin::Expression(2) {}
void eval(dolfin::Array<double>& values, const dolfin::Array<double>& coords) const {
// shift point of origin
double x = coords[0]-0.5;
double y = coords[1];
double z = coords[2]-0.5;
// transform to polar coordinates
double r = std::sqrt(std::pow(x,2) + std::pow(y,2) + std::pow(z,2));
// calculate distance to contourline at phi=0
double xi = r - 0.25;
// calculate phi
if (std::abs(xi) < M_PI/2.0*0.02) {
values[0] = std::sin(xi/0.02);
} else {
if (xi > 0.0) values[0] = +1.0;
if (xi < 0.0) values[0] = -1.0;
}
// set mu
values[1] = 0.0;
}
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<F, std::shared_ptr<F>, dolfin::Expression>
(m, "F")
.def(py::init<>());
}
"""
# create mesh
comm_osh = omega_h.world()
mesh_osh = omega_h.build_box(comm_osh, omega_h.Family.SIMPLEX, 1.0, 1.0, 1.0, 10, 10, 10)
mesh_osh.set_parting(omega_h.Parting.GHOSTED, 1)
omega_h.add_implied_metric_tag(mesh_osh)
mesh_osh.set_parting(omega_h.Parting.ELEM_BASED, 0)
mesh = dolfin.Mesh()
omega_h.mesh_to_dolfin(mesh, mesh_osh)
# create function
P1 = dolfin.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
CH = dolfin.FunctionSpace(mesh, dolfin.MixedElement(P1, P1))
V = dolfin.Function(CH)
v_init = dolfin.CompiledExpression(dolfin.compile_cpp_code(f_code).F(), degree=1)
V.interpolate(v_init)
# output
file = dolfin.XDMFFile(dolfin.MPI.comm_world, "output/phi_3d.xdmf")
file.parameters["rewrite_function_mesh"] = True
file.parameters["flush_output"] = True
phi, _ = V.split()
phi.rename(PHASEFIELD, PHASEFIELD)
file.write(phi, 0)
i = 0
while i<10:
i+=1
dolfin.info("step: {}".format(i))
# adapt
omega_h.function_from_dolfin(mesh_osh, phi._cpp_object, PHASEFIELD)
mesh_osh.set_parting(omega_h.GHOSTED, 1)
metric_input = omega_h.MetricInput()
source = omega_h.MetricSource(omega_h.VARIATION, 1e-2, PHASEFIELD)
metric_input.add_source(source)
metric_input.should_limit_lengths = True
metric_input.max_length = 2
metric_input.min_length = 0
metric_input.should_limit_gradation = True
metric_input.should_limit_element_count = True
metric_input.max_element_count = 5e5
omega_h.generate_target_metric_tag(mesh_osh, metric_input)
opts = omega_h.AdaptOpts(mesh_osh)
opts.min_quality_allowed = 0.1
#opts.min_quality_desired = 0.2
opts.verbosity = omega_h.EACH_ADAPT
while omega_h.approach_metric(mesh_osh, opts):
omega_h.adapt(mesh_osh, opts)
# refined mesh
mesh_fine = dolfin.Mesh()
omega_h.mesh_to_dolfin(mesh_fine, mesh_osh)
P1 = dolfin.FiniteElement("Lagrange", mesh_fine.ufl_cell(), 1)
CH = dolfin.FunctionSpace(mesh_fine, dolfin.MixedElement(P1, P1))
V = dolfin.Function(CH)
v_init = dolfin.CompiledExpression(dolfin.compile_cpp_code(f_code).F(), degree=1)
V.interpolate(v_init)
phi, _ = V.split()
phi.rename(PHASEFIELD, PHASEFIELD)
file.write(phi, i)
I suspect the phi field has become very non-smooth in some part of the mesh... it would be interesting to look more closely at it. Typically gradation becomes slow when there are a few points in the mesh that request very very fine resolution due to spikes or noise in the input field. Another thing to experiment with is should_limit_element_count = False
, as that process also iterates and could be iterating long for some unforeseen reason. I am also concerned by the decreasing min_quality_allowed
. Omega_h (and your simulation) will perform much better if you can raise that 0.3
.
Thx @ibaned! I am still working on this issue and will leave this question open for future discussions if it is okay for you.