# Negative diffusion viscosity detected,bug about velocity boundary

**URL:** <https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620>\
**Category:** ASPECT\
**Created:** [September 9, 2026, 3:01am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620 "2026-09-09T03:01:05Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![guit-hub](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/guit-hub/32/3105_2.png) [@guit-hub](https://community.geodynamics.org/u/guit-hub)\
**Post date:** [September 9, 2026, 3:01am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/1 "2026-09-09T03:01:05Z")

</div>

Dear developers:

I hope you’re doing everything well.When I run a model like this:

 ![IMG_3373.PNG](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/d/d496d67ba8744e7194508d1550b8e723c99c3829.jpeg)  
I want set velocity boundary conditions like this:  
0-140km,left,zero velocity,right,v\_x=influx,v\_y=0  
140-450km,left and right,v\_y=0,v\_x keeps variable outflux to account for mass conservation.  
So I adopt geometry model:Box with lithosphere indicators:  
subsection Geometry model

set Model name = box with lithosphere boundary indicators

subsection Box with lithosphere boundary indicators

```
set Box origin X coordinate =0

set Box origin Y coordinate = 0

set Lithospheric thickness = 140e3

set Use merged grids = true

set X extent = 800e3 # 500 km

set Y extent = 450e3 # 450 km

set X repetitions = 16

set Y repetitions = 10

```

end

end  
And I set Boundary velocity model like this:  
subsection Boundary velocity model

set Prescribed velocity boundary indicators = \

```
left lithosphere:zero velocity,\\

right lithosphere:function, \\

left y:function, \\

right y:function

```

set Tangential velocity boundary indicators = bottom

# set Zero velocity boundary indicators =

subsection Function

```
set Coordinate system = cartesian

set Function constants = influx=-10e-3 #10mm/yr

set Variable names = x,y

set Function expression = influx;0

```

end

end  
And correspondingly Temperature boundary:  
subsection Boundary temperature model

set Fixed temperature boundary indicators = top, bottom

set List of model names = box with lithosphere boundary indicators

subsection Box with lithosphere boundary indicators

```
\# set Boundary indicator to temperature mappings = top:300,bottom:1610

set Bottom temperature = 1610

set Top temperature = 300

```

end

end  
Then the “Negative diffusion viscosity detected” come up with the first timestep when steping to solve Stokes system  
`An error occurred in line <124> of file </hpcdata/home/research/hzhang02/opt/aspect-3.0.0/source/material_model/rheology/diffusion_creep.cc> in function `  
` double aspect::MaterialModel::Rheology::DiffusionCreep<dim>::compute_viscosity(double, double, double, unsigned int, const std::vector<double>&, const std::vector<unsigned int>&) const [with int dim = 2]`  
`The violated condition was: `  
` viscosity_diffusion > 0.0`  
`Additional information: `  
` Negative diffusion viscosity detected. This is unphysical and should`  
` not happen. Check for negative parameters. Temperature and pressure`  
` are 300 K, -431410827152.85846 Pa.`  
I think Temperature has no problem because it’s just simple and the temperature in the bug information is 300K,but the pressure is singular.And the material model is also normal because I can run it smoothing without influx boundary.  
So I think there’s something wrong with Boundary model.Could you offer me some advice? Any guidance is usefull for me the fresh.Thanks!  
Then I just try to prescribe the velocity boundary mor precisely,like this:  
subsection Boundary velocity model

set Prescribed velocity boundary indicators = \

```
left lithosphere:zero velocity,\\

right lithosphere x:function, \\

right lithosphere y: zero velocity, \\

left y:zero velocity, \\

right y:zero velocity

```

set Tangential velocity boundary indicators = bottom

# set Zero velocity boundary indicators =

subsection Function

```
set Coordinate system = cartesian

set Function constants = influx=-10e-3 #10mm/yr

set Variable names = x,y

set Function expression = influx

```

end

end  
Then some grammer bug happens:  
`An error occurred in line <340> of file </hpcdata/home/research/hzhang02/opt/aspect-3.0.0/source/boundary_velocity/interface.cc> in function`  
`void aspect::BoundaryVelocity::Manager::parse_parameters(dealii::ParameterHandler&) [with int dim = 2]`  
`The violated condition was:`  
`boundary_velocity_indicators[boundary_id].first == comp`  
`Additional information:`  
`Different velocity plugins for the same boundary have to have the same`  
`component selector. This was not the case for boundary: right`  
`lithosphere, for plugin: zero velocity, with component selector: y`

---

<div class="post-metadata">

**Author:** ![guit-hub](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/guit-hub/32/3105_2.png) [@guit-hub](https://community.geodynamics.org/u/guit-hub)\
**Post date:** [September 9, 2026, 3:35am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/2 "2026-09-09T03:35:37Z")

</div>

Best wishes,

Weidong Pian

---

<div class="post-metadata">

**Author:** ![daniel.douglas](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/daniel.douglas/32/2126_2.png) [@daniel.douglas](https://community.geodynamics.org/u/daniel.douglas)\
**Post date:** [September 9, 2026, 11:35am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/3 "2026-09-09T11:35:33Z")

</div>

Hi Weidong,

Thanks for posting on the forum and for showing that helpful figure! I agree that it likely isn’t the temperature boundary conditions, and the good news is that your stated goal:

> [@guit-hub](#):
>
> want set velocity boundary conditions like this:  
> 0-140km,left,zero velocity,right,v\_x=influx,v\_y=0  
> 140-450km,left and right,v\_y=0,v\_x keeps variable outflux to account for mass conservation.  
> So I adopt geometry model:Box with lithosphere indicators:

seems to me like it should be possible, and also that you are already using the geometry model (Box with lithosphere boundary indicators) that should allow you to do this.

The formatting of your copy/pasted text ended up a little wonky so I might be missing something, but from what I can tell this part of your your prm file:

```auto
subsection Boundary velocity model

set Prescribed velocity boundary indicators = left lithosphere:zero velocity,\\
                                                                        right lithosphere:function, \\
                                                                        left y:function, \\
                                                                        right y:function

  subsection Function
    set Coordinate system = cartesian
    set Function constants = influx=-10e-3
    set Variable names = x,y
    set Function expression = influx;0
  end

  set Tangential velocity boundary indicators = bottom

end

```

looks correct to me. Two things that seem like they are missing are:

1. What is the boundary condition on the top boundary?
2. You can prescribe what depth–pressure profile used for the right/left x boundary condition in the [Boundary traction model](https://aspect-documentation.readthedocs.io/en/latest/parameters/Boundary_20traction_20model.html) by doing something like:

```auto
subsection Boundary traction model
  Prescribed traction boundary indicators = right x: initial lithostatic pressure, left x: initial lithostatic pressure
  subsection Initial lithostatic pressure
  ......
  ......
  end
end

```

this might constrain the inflow/outflow more on these boundaries.  
3. Using a `Tangential velocity boundary` on the bottom, with open right/left x boundary conditions can introduce a horizontal nullspace that, if you aren’t removing, could cause unphysical mode behavior.

If you aren’t remove the nullspace already, take a look at the [Nullspace removal manual entry](https://aspect-documentation.readthedocs.io/en/latest/parameters/Nullspace_20removal.html#parameters-nullspace-20removal-remove-20nullspace) for what the syntax looks like. I would suspect that you want to remove the linear x momentum nullspace.

If this doesn’t help, for next steps it would be helpful to see your entire .prm if you wouldn’t mind uploading it.

Cheers,  
Daniel

---

<div class="post-metadata">

**Author:** ![guit-hub](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/guit-hub/32/3105_2.png) [@guit-hub](https://community.geodynamics.org/u/guit-hub)\
**Post date:** [September 10, 2026, 12:00am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/4 "2026-09-10T00:00:47Z")

</div>

Hi professor Daniel:

Thanks for your attention.Your advice just enlight some key points that I’ve never considered. The top boundary is set as free surface because I want to see surface topography. I have uploaded my parameters file.Following your advice,I add removal of linear x momentum.

```cpp
subsection Nullspace removal

  set Remove nullspace = linear x momentum

end

```

But it just still report the same problem.Then I add lithostastic pressure:

```auto
subsection Boundary traction model

  set Prescribed traction boundary indicators = right x:initial lithostatic pressure,left x:initial lithostatic pressure

  #the start point of pressure profile

  subsection Initial lithostatic pressure

    set Number of integration points = 5000

    set Representative point = 0,0

  end

end

```

Finally,it runs up smoothly.  
But it’s still can’t be sure that everything is ok,because the velocity statistics seems to large

```auto
*** Timestep 17: t=3473.46 years, dt=193.828 years
Solving mesh displacement system... 5 iterations.
Solving temperature system... 13 iterations.
Solving upper_crust system ... 13 iterations.
Solving lower_crust system ... 14 iterations.
Solving Mantle_lithosphere system ... 15 iterations.
Solving asthenosphere system ... 13 iterations.
Solving weak_zone system ... 15 iterations.
Solving Stokes system (GMG)... 388+0 iterations.

Postprocessing:
Temperature min/avg/max: 300 K, 1451 K, 1610 K
RMS, max velocity: 0.437 m/year, 7.63 m/year
Pressure min/avg/max: -7.199e+07 Pa, 6.833e+09 Pa, 1.384e+10 Pa
Heat fluxes through boundary parts: -8.099e+05 W, -1.121e+06 W, -1107 W, 7.095e+05 W, -7.619e+04 W, -1.206e+05 W
Topography min/max: 0 m, 3864 m
Writing heat flux map output-hangai6/heat_flux.00017
Pressure at top/bottom of domain: 2.673e+07 Pa, 1.384e+10 Pa
Computing dynamic topography

```

Although the velocity goes down with time stepping.  
I’ll see the results after the job is done.

Thanks again,

Weidong Pian

[hangai1.prm](https://community.geodynamics.org/uploads/short-url/qmSzkXyuOHwBFloQHCmaFu5Wbn8.prm) (9.8 KB)

---

<div class="post-metadata">

**Author:** ![daniel.douglas](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/daniel.douglas/32/2126_2.png) [@daniel.douglas](https://community.geodynamics.org/u/daniel.douglas)\
**Post date:** [September 10, 2026, 6:30pm UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/5 "2026-09-10T18:30:42Z")

</div>

Hi Weidong,

Great to hear that the model is running now! Because you are using a free surface, the initial large velocities that decay with time do not necessarily mean that the model is wrong.  
Looking at the .prm you uploaded, you have an initial composition that specifies a lithospheric root in a subset of the model. This means that relative to the parts of your model with a lithosphere of “normal thickness”, the integrated buoyancy force of the lithosphere with a root will be different. Combined with the free surface, this can result in very large velocities initially while the model tries to isostatically compensate the lithospheric root.

As long as the velocities eventually decay to reasonable values (over say 10 - 100 kyr) then this is likely just due to initial isostatic adjustment of the free surface. Keep us posted on whether the results look right! If you are having trouble with solver behaviour during this adjustment period, it could be helpful to reduce the CFL number from 1.0 down to something like 0.1.

Cheers,  
Daniel

---

<div class="post-metadata">

**Author:** ![guit-hub](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/guit-hub/32/3105_2.png) [@guit-hub](https://community.geodynamics.org/u/guit-hub)\
**Post date:** [September 12, 2026, 2:36am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/6 "2026-09-12T02:36:52Z")

</div>

Hi Professor Douglas:

I’ve got my results and I got some questitions.

Firstly, When the model runs to ~6Ma,the velocity field just lose its control.And consequently,the temperature field just became strange too.

```auto
*** Timestep 962: t=6.92663e+06 years, dt=71.1765 years
Solving mesh displacement system... 5 iterations.
Solving temperature system... 14 iterations.
Solving upper_crust system ... 15 iterations.
Solving lower_crust system ... 16 iterations.
Solving Mantle_lithosphere system ... 17 iterations.
Solving asthenosphere system ... 17 iterations.
Solving weak_zone system ... 17 iterations.
Solving Stokes system (GMG)... 1000+37 iterations.

Postprocessing:
Temperature min/avg/max: 24.52 K, 1413 K, 6262 K
RMS, max velocity: 0.896 m/year, 18.3 m/year
Pressure min/avg/max: -1.133e+08 Pa, 6.837e+09 Pa, 1.398e+10 Pa
Heat fluxes through boundary parts: 3.252e+06 W, -5.112e+06 W, -2115 W, 3.685e+05 W, -1.939e+06 W, 1.631e+06 W
Topography min/max: -117.1 m, 3661 m
Writing heat flux map influx1/heat_flux.00962
Pressure at top/bottom of domain: -1.536e+07 Pa, 1.385e+10 Pa

```

I think these all result from the velocity discontinue point between the lithosphere sides and asthenosphere sides,which can be showed on the following snap of `Time=~6.9Ma`

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/d/d135df967fb26cd0b6fe9958827e8ba008238bb3.jpeg)

Finally the model ended up with Negative diffusion viscosity deteched again.

Maybe this problem dosen’t bother me,because I think it will not interfere the question I want to study.(I want reproduce the performance reported in (Comeau et al., 2021),which shows how convergence and phase transition imply on lithosphere removal. **But I do want to know hwo to avoid this discontinuity out of curiosity.**

Secondly,I tried control some parameters including different influx,a shallower phase transition depth(35km) to run different model,but **I failed to see that materal of lower crust sinks into asthenosphere after phase transition for a denser density.**

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/2/2372fc691c413da87148ca4b04a05cfec85bb467.png)

So I refer to Compositing model to get a more “complex” model,like this:

```auto
# ------------------------------------------------------------

# Material model: Visco-plastic with phase transition

# ------------------------------------------------------------

subsection Material model

  # set Model name = visco plastic

  set Model name = compositing

  set Material averaging = arithmetic average

  subsection Visco Plastic

    # Viscosity parameters (dislocation creep)

    set Prefactors for dislocation creep = 1.76e-19, 0.876e-27, 0.556e-18, 2.19e-15, 1.76e-14, 2e-10

    set Stress exponents for dislocation creep = 1,4,4.7,3,3, 4

    set Activation energies for dislocation creep = 0,223e3, 485e3, 430e3, 430e3, 223e3 # J/mol

    set Activation volumes for dislocation creep = 0, 0, 0, 15e-6, 15e-6, 0 # m^3/mol

    set Reference strain rate = 1e-15 # 1/s

    set Reference temperature = 1600

    set Minimum viscosity = 1e17

    set Maximum viscosity = 1e24

    # Yield stress (Mohr-Coulomb)

    set Cohesions = 0,20e6, 0, 0, 0, 20e6 # Pa

    set Angles of internal friction = 12, 12, 12, 12, 12, 12 # degrees

    set Viscosity averaging scheme = arithmetic

    # Phase transition (eclogitization)

    # set Number of phases = 1

    set Phase transition depths = 35e3 # m

    set Phase transition temperatures = 1110 # K

    set Phase transition Clapeyron slopes = 1e6 # Pa/K (1 MPa/K)

    set Phase transition widths = 10e3 # m

    # set Phase transition density jumps = 400 # kg/m^3 (density increase)

    # set Corresponding phase for density jump = 0

    # Compositional dependencies

    # set Composition viscosity prefactor = 0

    # set Composition thermal expansivity = 0

    # set Composition thermal conductivity = 0

    set Thermal expansivities = 3e-5

    set Thermal conductivities = 0

    set Thermal diffusivities = 1.5e-5, 1.1e-6, 1.0e-6,5.5e-7,1.5e-5, 1e-6

    set Heat capacities = 1250, 750, 750, 1250, 1250, 750 # J/kg/K

    set Densities = 3300, 2700, 2800, 3300, 3300, 2800 # kg/m^3

  end

  subsection Latent heat

    set Phase transition Clapeyron slopes = 1e6 # Pa/K (1 MPa/K)

    set Phase transition density jumps = 400 # kg/m^3 (density increase)

    set Phase transition depths = 35e3 # m

    set Phase transition temperatures = 1110 # K

    set Phase transition widths = 10e3 # m

    set Reference density = 3250

    set Corresponding phase for density jump = 2

  end

  subsection Compositing

    set Density = visco plastic

    set Reaction terms = latent heat

    set Specific heat = visco plastic

    set Thermal conductivity = visco plastic

    set Thermal expansion coefficient = visco plastic

    set Viscosity = visco plastic

  end

end

```

I choose to do so,just because I want to avoid write a plugin.And `cookbooks/mantle_convection_with_continents_in_annulus/modelR.prm`is my reference.

```auto
An error occurred in line <773> of file </hpcdata/home/research/hzhang02/opt/aspect-3.0.0/include/aspect/plugins.h> in function
    static std::unique_ptr<_Tp> aspect::internal::Plugins::PluginList<InterfaceClass>::create_plugin(const string&, const string&) [with InterfaceClass = aspect::MaterialModel::Interface<2>; std::string = std::__cxx11::basic_string<char>]
The violated condition was: 
    name != "unspecified"
Additional information: 
    A plugin must have a name!
    
    This function was asked to create a plugin but no name for the plugin
    was provided. This may be due to the fact that you did not explicitly
    specify a name for this plugin in your input file and ASPECT does not
    provide a default for this kind of plugin, for example because no
    generally useful plugin exists. An example is that there is no default
    geometry: You need to explicitly provide one in the input file, and it
    seems like you have not done so.
    
    To find out which kind of plugin this function tries to create, take a
    look at the backtrace of this error message.
    
    The place that called this function also provided as additional
    information this:
    
    <Material model::Model name>

Stacktrace:
-----------
#0 aspect: aspect::internal::Plugins::PluginList<aspect::MaterialModel::Interface<2> >::create_plugin(std:: __cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__ cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&)
#1 aspect: std::unique_ptr<aspect::MaterialModel::Interface<2>, std::default_delete<aspect::MaterialModel::Interface<2> > > aspect::MaterialModel::create_material_model<2>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&)
#2 aspect: aspect::MaterialModel::Compositing<2>::parse_parameters(dealii::ParameterHandler&)
#3 aspect: aspect::Simulator<2>::Simulator(ompi_communicator_t*, dealii::ParameterHandler&)
#4 aspect: void run_simulator<2>(std:: __cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, std::__ cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, bool, bool, bool, bool)
#5 aspect: main
--------------------------------------------------------

Aborting!

```

I have no idea:  
**1.Why ASPECT reports this error?Do I really have to write a plugin?It seems like the origin modelR.prm dosen’t have this problem.If I have to do so,I’ll take time to make it.  
2. If I have set relative parameters about phase transition and other parameters in subsection Latent heat,do I need set those above again in subsection visco plastic?If I just set those value in both subsection like what it dose in modelR.prm,which one the ASPECT will choose?**

---

<div class="post-metadata">

**Author:** ![daniel.douglas](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/daniel.douglas/32/2126_2.png) [@daniel.douglas](https://community.geodynamics.org/u/daniel.douglas)\
**Post date:** [September 14, 2026, 4:26pm UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/7 "2026-09-14T16:26:59Z")

</div>

Hi Weidong,

Great to see that the model is running for millions of years! The crash at just under 7 Myr is tricky, one thing that might be worth trying is discontinuous discretization? That might help with the change in the boundary conditions across the lithosphere/asthenosphere boundary.

For your first question, you should not have to create your own plugin and I suspect that this error is being shown because you are missing some of the required entries of the `Compositing` material model. Specifically, you are missing `Entropy derivative temperature` and `Entropy derivative pressure`, which are defined in the modelR.prm example you reference.

As you point out, the Visco Plastic material model can deal with phase transitions without being composited with another material model, so unless you need the latent heat material model it would be easier to just use Visco Plastic.

To answer your seceond question, the `Compositing` material model will use the values from the base material model that you specify. So in the example you posted, the density will come from `visco plastic`, even if you define the density in `latent heat`. The reaction terms will come from `latent heat`, etc.,

Cheers,  
Daniel

---

<div class="post-metadata">

**Author:** ![guit-hub](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/guit-hub/32/3105_2.png) [@guit-hub](https://community.geodynamics.org/u/guit-hub)\
**Post date:** [September 15, 2026, 3:08am UTC](https://community.geodynamics.org/t/negative-diffusion-viscosity-detected-bug-about-velocity-boundary/4620/8 "2026-09-15T03:08:28Z")

</div>

Hi Professor Douglas:

Thanks for your kindly enlightment to me.Now the model works again,I’ll see what happens next.I’ve got some basic understanding of `Visco plastic model`from many papers.Also it’s nessary to learn about `Latent model` to deduce whether I have to use that model,because I try not to make my model too complex.

Thanks again!The kindly atmosphere of the community alway keeps one studying and searching.

Best wishes,

Weidong Pian