Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

warning: gradation limiting is up to step 2200 #246

Open
bonh opened this issue Aug 1, 2018 · 2 comments
Open

warning: gradation limiting is up to step 2200 #246

bonh opened this issue Aug 1, 2018 · 2 comments
Labels

Comments

@bonh
Copy link

bonh commented Aug 1, 2018

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:

  1. Start with a very coarse and uniform mesh (mesh1)
  2. Interpolate the initial conditions on mesh1
  3. Adapt mesh1 and get mesh2
  4. Interpolate the initial conditions on mesh2
  5. Repeat steps 3 and 4 a few times
  6. 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)
@ibaned
Copy link
Collaborator

ibaned commented Aug 1, 2018

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.

@bonh
Copy link
Author

bonh commented Sep 3, 2018

Thx @ibaned! I am still working on this issue and will leave this question open for future discussions if it is okay for you.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants