Coding dirichletbc from material property #23797

vermaprk · 2023-03-21T06:47:31Z

vermaprk
Mar 21, 2023

I have written the include and src files for executing Dirichlet boundary condition taking values from material property. This runs well on simple problems. However, this is leading to convergence issues for complex problems. I have used IntegratedBC as the parent class. Please look into my code (attached)

header file

#pragma once

#include "IntegratedBC.h"

class EgMatDirichlet : public IntegratedBC
{
public:
  EgMatDirichlet(const InputParameters & parameters);

  static InputParameters validParams();

protected:
  virtual Real computeQpResidual() override;

private:
  const MaterialProperty<Real> & _boundary_prop;
};

src file

#include "EgMatDirichlet.h"

registerMooseObject("goldfishApp", EgMatDirichlet);

InputParameters
EgMatDirichlet::validParams()
{
  InputParameters params = IntegratedBC::validParams();

  // Specify input parameters that we want users to be able to set:
  params.addRequiredParam<MaterialPropertyName>("boundary_material", "Value on the boundary");
  return params;
}

EgMatDirichlet::EgMatDirichlet(const InputParameters & parameters)
  : IntegratedBC(parameters),
    // store the user-specified parameters from the input file...
    _boundary_prop(getMaterialProperty<Real>("boundary_material"))
    
{
}

Real
EgMatDirichlet::computeQpResidual()
{
  
  return _boundary_prop[_qp] - _u[_qp];
}

Answered by vermaprk

Mar 23, 2023

Aren't both same, as the residue has to be near to zero.
I tried the following change on simplediffusion.i.

ADReal
ADEgMatDirichlet::computeQpResidual()
{
  
  return _u[_qp] - _boundary_prop[_qp];
}

The postprocessor gives:
Postprocessor Values:

+----------------+----------------+----------------+
| time           | left_c         | right_c        |
+----------------+----------------+----------------+
:                :                :                :
|   1.785000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.795000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.815000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.835000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.855000e+02 |   1.00…

View full answer

GiudGiud · 2023-03-21T15:51:09Z

GiudGiud
Mar 21, 2023
Collaborator

hello

You are missing the contribution to the Jacobian. Does it converge with PJFNK for the solve_type?
You should just use AD to make this boundary condition and it will handle that part if you want to make it work with Newton

Guillaume

6 replies

maxnezdyur Mar 22, 2023

It looks like you are using an integrated BC class to set a Dirichlet BC. I think you need a penalty parameter to have the values match closer. You will also need to multiply by the test function. Try the below:
return _penalty * _test[_i][_qp] * (_u[_qp] - _boundary_prop[_qp]); where _penalty is a large number like 1e7 (problem dependent).

vermaprk Mar 22, 2023
Author

Yes, I am using ADIntegratedBC class as for ADNodalBC I couldn't find the getADMaterialProperty<Real> as a Public Member Function.
Is there a way to do this for ADNodalBC ? (I have a limited knowledge of C++)
I tried your method, but this also leads to non-convergence for the simple diffusion problem

ADReal
ADEgMatDirichlet::computeQpResidual()
{
  Real penalty = 1e7;
  return penalty * _test[_i][_qp] * (_u[_qp] - _boundary_prop[_qp]);
}

my input file for simple diffusion:

[Mesh]
    type = GeneratedMesh
    dim = 1
    xmin = 0
    xmax = 200e-6
    nx = 200
[]
  
[Variables]
    [c]
    []
[]

# Initial conditions.
[ICs]
    [conc]
        variable = c
        type = FunctionIC
        function = ic_func_c
    []
[]

[BCs]
    [left_c]
        type = DirichletBC
        variable = c
        boundary = 'left'
        value = 1000
    []
    [right_c]
        type = ADEgMatDirichlet
        variable = c
        boundary = 'right'
        boundary_material=bndconst
    []
[]
  
[Functions]
    [ic_func_c]
        type = ParsedFunction
        expression = 1000
    []
[]
  
[Kernels]
    [cdot]
        type = TimeDerivative
        variable = c
    []
    [c_diff]
        type=MatDiffusion
        diffusivity= 7.5e-11
        variable = c
    []
[]


[Materials]
    [matconstant]
        type = ADGenericConstantMaterial
        prop_names = 'bndconst'
        prop_values = '300'
    []
[]

[Executioner]
    type = Transient
    solve_type = 'PJFNK'
    dt=2
    end_time = 200
[]

[Postprocessors]
    [left_c]
      type = PointValue
      point = '0 0 0'
      variable = c
    []
    [right_phi]
      type = PointValue
      point = '200e-6 0 0'
      variable = c
    []
  
  []
[Outputs]
    exodus = true
[]

GiudGiud Mar 22, 2023
Collaborator

Material properties aren’t well defined on nodes because nodes can be shared between subdomains and we wouldn’t know which subdomain to pick to evaluate the property (their expression can vary by subdomain)

hence why there is not getMatProp on nodal classes

now that you are using AD, does Newton converge with the old integrated BC?

vermaprk Mar 23, 2023
Author

Thanks for the help.
Newton converges for the diffusion problem. However the constant material property value 300 is showing in postprocessor as 2.999999e+02 :

time           | left_c         | right_c        |
+----------------+----------------+----------------+
:                :                :                :
|   1.905000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.907500e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.912500e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.915000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.920000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.925000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.935000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.945000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.950000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.960000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.980000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.990000e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.992500e+02 |   1.000000e+03 |   2.999999e+02 |
|   1.995000e+02 |   1.000000e+03 |   2.999999e+02 |
|   2.000000e+02 |   1.000000e+03 |   2.999999e+02 |

However, in case of my original problem newton doesn't converges at all.
PJFNK works well with DirichletBC and the computation time is almost double for ADEgMatDirichlet

GiudGiud Mar 23, 2023
Collaborator

in ADEgMatDirichlet, how do you retrieve the material property?

Is this an AD material property? If not you re still missing terms in the Jacobian

vermaprk · 2023-03-23T05:06:53Z

vermaprk
Mar 23, 2023
Author

ADEgMatDirichlet.C

#include "ADEgMatDirichlet.h"

registerMooseObject("goldfishApp", ADEgMatDirichlet);

InputParameters
ADEgMatDirichlet::validParams()
{
  InputParameters params = ADIntegratedBC::validParams();

  // Specify input parameters that we want users to be able to set:
  params.addRequiredParam<MaterialPropertyName>("boundary_material", "Value on the boundary");
  return params;
}

ADEgMatDirichlet::ADEgMatDirichlet(const InputParameters & parameters)
  : ADIntegratedBC(parameters),
    // store the user-specified parameters from the input file...
    _boundary_prop(getADMaterialProperty<Real>("boundary_material"))
    
{
}

ADReal
ADEgMatDirichlet::computeQpResidual()
{
  
  return _boundary_prop[_qp] - _u[_qp];
}

And the input file material block is:

[Materials]
    [matconstant]
        type = ADGenericConstantMaterial
        prop_names = 'bndconst'
        prop_values = '300'
    []
[]

1 reply

GiudGiud Mar 23, 2023
Collaborator

shouldn't it be
_u[_qp] - _boundary_prop[_qp];

vermaprk · 2023-03-23T08:57:35Z

vermaprk
Mar 23, 2023
Author

Aren't both same, as the residue has to be near to zero.
I tried the following change on simplediffusion.i.

ADReal
ADEgMatDirichlet::computeQpResidual()
{
  
  return _u[_qp] - _boundary_prop[_qp];
}

The postprocessor gives:
Postprocessor Values:

+----------------+----------------+----------------+
| time           | left_c         | right_c        |
+----------------+----------------+----------------+
:                :                :                :
|   1.785000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.795000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.815000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.835000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.855000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.875000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.895000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.915000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.925000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.945000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.965000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.975000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.985000e+02 |   1.000000e+03 |   3.000001e+02 |
|   1.992500e+02 |   1.000000e+03 |   3.000001e+02 |
|   2.000000e+02 |   1.000000e+03 |   3.000001e+02 |
+----------------+----------------+----------------+

Interestingly, this change made in ADEgMatDirichlet leads to convergence issue on the original problem at the end time steps for PJFNK. Before, this change the code was running completely with PJFNK.

Solve Did NOT Converge!
Aborting as solve did not converge

Time Step 177, time = 199.999, dt = 1e-13
 0 Nonlinear |R| = 4.105825e-06
      0 Linear |R| = 4.105825e-06
      1 Linear |R| = 8.327757e-13
Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH iterations 0
 Solve Did NOT Converge!
Aborting as solve did not converge


*** ERROR ***
The following error occurred in the object "TimeStepper", of type "ConstantDT".

Solve failed and timestep already at or below dtmin, cannot continue!

3 replies

GiudGiud Mar 23, 2023
Collaborator

This is not a real divergence.

Time Step 177, time = 199.999, dt = 1e-13
0 Nonlinear |R| = 4.105825e-06
0 Linear |R| = 4.105825e-06
1 Linear |R| = 8.327757e-13

you're converged, but the solver tries to take a 1e-13 timestep and it throws everything off

if you increase the timestep_tol to 1e-12 this will be fixed

Could you please paste the log for the Newton solve_type?

vermaprk Mar 24, 2023
Author

Please find the log file for newton attached.
nohup.txt

input file:

alpha = 0.5
cinit = 1000
cmax = 3920
DLi = 7.5e-11
R_const = 8.314
F_const = 96487
ka = 1e-10
t0 = 0.36
I_app = -10
T = 300

[Mesh]
    type = GeneratedMesh
    dim = 1
    xmin = 0
    xmax = 200e-6
    nx = 200
[]

[Variables]
    [c]
    []
    [phi]
    []
[]

# Initial conditions.
[ICs]
    [conc]
        variable = c
        type = FunctionIC
        function = ic_func_c
    []
[]

[BCs]
    [left_c]
        type = ADMatNeumannBC
        variable = c
        boundary = 'left'
        value = 1
        boundary_material = Nflux
    []
    # [right_c]
    #     type = ADMatNeumannBC
    #     variable = c
    #     boundary = 'right'
    #     value = -1
    #     boundary_material = Nflux
    # []
    [right_c]
        type = ADEgMatDirichlet
        variable = c
        boundary = 'right'
        boundary_material=bndconst
    []
    # [right_c]
    #         type = ADDirichletBC
    #         variable = c
    #         boundary = 'right'
    #         value = 300
    #     []
    [left_phi]
        boundary = 'left'
        type = MatDirichletBC
        variable = phi
        alpha = '${alpha}'
        Iapp = '${I_app}'
        cmax = '${cmax}'
        temp = '${T}'
        ka = '${ka}'
        some_var = c
    []
[]

[AuxVariables]
    [velocity]
        family = LAGRANGE_VEC
    []
[]

[AuxKernels]
    [vec]
        type = VectorFunctionAux
        variable = velocity
        function = velocity_func
        execute_on = 'INITIAL TIMESTEP_END'
    []
[]

[Functions]
    [velocity_func]
        type = ParsedVectorFunction
        expression_x = '1'
    []

    [ic_func_c]
        type = ParsedFunction
        expression = '${cinit}'
    []
[]

[Kernels]
    [cdot]
        type = TimeDerivative
        variable = c
    []
    [c_diff]
        type=ADMatDiffusion
        diffusivity='${DLi}'
        variable = c
    []
    [c_mig]
        type=ExampleMigration
        variable = phi
        velocity = velocity
        reaction_rate = coeff_mig
        v=c
    []
    [advection]
        type = LevelSetAdvection
        velocity = velocity
        variable = phi
    []

    [rhs]
        type = ADMaterialPropertyValue
        variable = phi
        prop_name = IbyK
    []
[]

[Materials]
    [matconstant]
        type = ADGenericConstantMaterial
        prop_names = 'bndconst'
        prop_values = '1000'
    []
    [Boundaryflux]
        type = ADParsedMaterial
        property_name = Nflux
        expression = '${I_app}/${F_const}'
    []
    [migrationcoeff]
        type = ADParsedMaterial
        property_name = coeff_mig
        expression = '(2*${R_const}*${T}*(1-${t0}))/(${F_const})'
    []
    [ionicconductivity]
        type = ADParsedMaterial
        coupled_variables = 'c'
        property_name = kappa
        expression = '0.1147*exp(520/${R_const}*(1/298-1/${T}))*(c/1000)^3
                        -2.238*exp(1010/${R_const}*(1/298-1/${T}))*(c/1000)^1.5
                        +2.915*exp(1270/${R_const}*(1/298-1/${T}))*(c/1000)'
        outputs = exodus
    []
    [IbyK]
        type = ADParsedMaterial
        property_name = IbyK
        material_property_names = 'kappa'
        expression = '${I_app}/kappa'
        outputs = exodus
    []
[]

[Executioner]
    type = Transient
    solve_type = 'Newton'
    # solve_type = 'PJFNK'
    dt = 0.2
    end_time = 200
[]

[Postprocessors]
    [left_phi]
      type = PointValue
      point = '0 0 0'
      variable = phi
    []
    [right_phi]
      type = PointValue
      point = '200e-6 0 0'
      variable = phi
    []
    [right_c]
        type = PointValue
        point = '200e-6 0 0'
        variable = c
    []
    [I0]
        type = ParsedPostprocessor
        function = '${F_const}*${ka}*pow((${cmax}-right_c),${alpha})*pow(right_c,${alpha})'
        pp_names = 'right_c'
    []
    [voltage]
        type = ParsedPostprocessor
        function = 'right_phi + (${R_const}*${T})/(${alpha}*${F_const})*(asinh(${I_app}/I0))'
        pp_names = 'right_phi I0'
    []
[]


[Outputs]
    exodus = true
    [./csv]
        type = CSV
    [../]
[]

maxnezdyur Mar 24, 2023

Try the code and input file below.

I think you need to scale your problem. I think there are some conditioning problems.
I think that for an integrateBC to model a DirichletBC like you want it needs to be in the form written below
In my previous post, I said the penalty term is problem dependent. Using 1e2 as the penalty parameter, the solution converges and is accurate (post processor attached below).
I think the reason your problem worked somewhat for _u[_qp] - _boundary_prop[_qp] as the BC is because of the parameters you had, but as I increased the diffusivity, the BC below would still be accurate but the _u[_qp] - _boundary_prop[_qp] method would get progressively worse results.

#include "EgMatDirichlet.h"

registerMooseObject("goldfishApp", EgMatDirichlet);

InputParameters
EgMatDirichlet::validParams()
{
  InputParameters params = ADIntegratedBC::validParams();

  // Specify input parameters that we want users to be able to set:
  params.addRequiredParam<MaterialPropertyName>("boundary_material", "Value on the boundary");
  params.addRequiredParam<Real>("penalty", "Penalty to enforce the boundary condition.");
  return params;
}

EgMatDirichlet::EgMatDirichlet(const InputParameters & parameters)
  : ADIntegratedBC(parameters),
    // store the user-specified parameters from the input file...
    _boundary_prop(getADMaterialProperty<Real>("boundary_material")),
    _penalty(getParam<Real>("penalty"))
{
}

ADReal
EgMatDirichlet::computeQpResidual()
{
  return _penalty * _test[_i][_qp] * (_u[_qp] - _boundary_prop[_qp]);
}

#pragma once

#include "ADIntegratedBC.h"

class EgMatDirichlet : public ADIntegratedBC
{
public:
  EgMatDirichlet(const InputParameters & parameters);

  static InputParameters validParams();

protected:
  virtual ADReal computeQpResidual() override;

private:
  const ADMaterialProperty<Real> & _boundary_prop;
  Real _penalty;
};

[Mesh]
    type = GeneratedMesh
    dim = 1
    xmin = 0
    xmax = 200e-6
    nx = 200
[]

[Variables]
    [c]
    []
[]

# Initial conditions.
[ICs]
    [conc]
        variable = c
        type = FunctionIC
        function = ic_func_c
    []
[]

[BCs]
    [left_c]
        type = DirichletBC
        variable = c
        boundary = 'left'
        value = 1000
    []
    [right_c]
        type = EgMatDirichlet
        variable = c
        boundary = 'right'
        boundary_material=bndconst
        penalty = 1e2
    []
[]

[Functions]
    [ic_func_c]
        type = ParsedFunction
        expression = 1000
    []
[]

[Kernels]
    [cdot]
        type = TimeDerivative
        variable = c
    []
    [c_diff]
        type=MatDiffusion
        diffusivity= 7.5e-11
        variable = c
    []
[]


[Materials]
    [matconstant]
        type = ADGenericConstantMaterial
        prop_names = 'bndconst'
        prop_values = '300'
    []
[]

[Executioner]
    type = Transient
    solve_type = 'NEWTON'
    dt=2
    end_time = 200
    nl_abs_tol = 1e-10
[]

[Postprocessors]
    [left_c]
      type = PointValue
      point = '0 0 0'
      variable = c
    []
    [right_phi]
      type = PointValue
      point = '200e-6 0 0'
      variable = c
    []

  []
[Outputs]
    exodus = true
[]

Time Step 100, time = 200, dt = 2
 0 Nonlinear |R| = 2.257844e-06
      0 Linear |R| = 2.257844e-06
      1 Linear |R| = 5.617709e-20
 1 Nonlinear |R| = 1.234283e-12
 Solve Converged!

Postprocessor Values:
+----------------+----------------+----------------+
| time           | left_c         | right_phi      |
+----------------+----------------+----------------+
:                :                :                :
|   1.720000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.740000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.760000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.780000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.800000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.820000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.840000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.860000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.880000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.900000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.920000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.940000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.960000e+02 |   1.000000e+03 |   3.000000e+02 |
|   1.980000e+02 |   1.000000e+03 |   3.000000e+02 |
|   2.000000e+02 |   1.000000e+03 |   3.000000e+02 |
+----------------+----------------+----------------+

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Coding dirichletbc from material property #23797

{{title}}

Replies: 3 comments 10 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Coding dirichletbc from material property #23797

vermaprk Mar 21, 2023

Replies: 3 comments · 10 replies

GiudGiud Mar 21, 2023 Collaborator

maxnezdyur Mar 22, 2023

vermaprk Mar 22, 2023 Author

GiudGiud Mar 22, 2023 Collaborator

vermaprk Mar 23, 2023 Author

GiudGiud Mar 23, 2023 Collaborator

vermaprk Mar 23, 2023 Author

GiudGiud Mar 23, 2023 Collaborator

vermaprk Mar 23, 2023 Author

GiudGiud Mar 23, 2023 Collaborator

vermaprk Mar 24, 2023 Author

maxnezdyur Mar 24, 2023

vermaprk
Mar 21, 2023

Replies: 3 comments 10 replies

GiudGiud
Mar 21, 2023
Collaborator

vermaprk Mar 22, 2023
Author

GiudGiud Mar 22, 2023
Collaborator

vermaprk Mar 23, 2023
Author

GiudGiud Mar 23, 2023
Collaborator

vermaprk
Mar 23, 2023
Author

GiudGiud Mar 23, 2023
Collaborator

vermaprk
Mar 23, 2023
Author

GiudGiud Mar 23, 2023
Collaborator

vermaprk Mar 24, 2023
Author