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

Jacobian tester for local assemblers. #2238

Merged
merged 12 commits into from Oct 21, 2018

Conversation

Projects
None yet
4 participants
@chleh
Copy link
Collaborator

chleh commented Oct 10, 2018

As requested, here comes a simple implementation of the Jacobian tester.

Configuration in the prj file:

<OpenGeoSysProject>
    <processes>
        <process>
            <jacobian_assembler>
                <type>CompareJacobians</type>
                <jacobian_assembler>
                    <type>CentralDifferences</type>
                    <component_magnitudes>1 1</component_magnitudes>
                    <relative_epsilons>1e-8 1e-8</relative_epsilons>
                </jacobian_assembler>
                <reference_jacobian_assembler>
                    <type>CentralDifferences</type>
                    <component_magnitudes>1 1</component_magnitudes>
                    <relative_epsilons>1e-6 1e-6</relative_epsilons>
                </reference_jacobian_assembler>
                <abs_tol>1e-18</abs_tol>
                <rel_tol>1e-7</rel_tol>
                <fail_on_error>true</fail_on_error>
                <log_file>/tmp/test.log</log_file>
            </jacobian_assembler>
Sample (Python) log file (Click to expand)
#!/usr/bin/env python
import numpy as np
from numpy import nan


### counter: 1 (begin)
# absolute tolerance of 1e-18 exceeded and relative tolerance of 1e-07 exceeded

counter = 1
num_dof = 8
abs_tol = 1e-18
rel_tol = 1e-07

local_x = np.array([2, 1, 1, 2, 1, 0, 0, 1])
local_x_dot = np.array([10, 0, 0, 10, 10, 0, 0, 10])
dxdot_dx = 10
dx_dx = 1

Jacobian_1 = np.array([
    [1.93984028211111e-10, -2.78579862194445e-11, -8.88539932722222e-11, -5.28579862194445e-11, -1.66666679687575e-11, 1.66666686149923e-11, 8.33333317658516e-12, -8.33333511528972e-12],
    [-6.11913195527778e-11, 2.27317361544444e-10, -3.61913195527778e-11, -1.05520659938889e-10, -1.66666663531703e-11, 1.66666650607006e-11, 8.33333188411546e-12, -8.33333446905487e-12],
    [-1.05520659938889e-10, -3.61913195527778e-11, 2.27317361544444e-10, -6.11913195527778e-11, -8.33333446905487e-12, 8.33333188411546e-12, 1.666666764564e-11, -1.66666650607006e-11],
    [-5.28579862194445e-11, -8.88539932722222e-11, -2.78579862194445e-11, 1.93984028211111e-10, -8.33333511528972e-12, 8.33333285346774e-12, 1.66666669994052e-11, -1.66666679687575e-11],
    [0, 0, 0, 0, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06],
    [0, 0, 0, 0, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06]])
Jacobian_2 = np.array([
    [1.93984028211111e-10, -2.78579862194445e-11, -8.88539932722222e-11, -5.28579862194445e-11, -1.66666666601319e-11, 1.66666666569007e-11, 8.33333333168153e-12, -8.33333335429975e-12],
    [-6.11913195527778e-11, 2.27317361544444e-10, -3.61913195527778e-11, -1.05520659938889e-10, -1.66666666472072e-11, 1.66666666472072e-11, 8.33333334137505e-12, -8.33333333168153e-12],
    [-1.05520659938889e-10, -3.61913195527778e-11, 2.27317361544444e-10, -6.11913195527778e-11, -8.33333334460623e-12, 8.33333334137505e-12, 1.66666666601319e-11, -1.66666666730566e-11],
    [-5.28579862194445e-11, -8.88539932722222e-11, -2.78579862194445e-11, 1.93984028211111e-10, -8.33333335753092e-12, 8.33333332845035e-12, 1.6666666643976e-11, -1.66666666569007e-11],
    [0, 0, 0, 0, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06],
    [0, 0, 0, 0, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06]])

# Jacobian_2 - Jacobian_1
abs_diff = np.array([
    [0, 0, 0, 0, 1.30862557845303e-18, -1.95809160627787e-18, 1.55096364853693e-19, 1.76098997594297e-18],
    [0, 0, 0, 0, -2.94036858368459e-19, 1.58650656548256e-18, 1.45725959477115e-18, 1.13737334226041e-18],
    [0, 0, 0, 0, 1.12444864518927e-18, 1.45725959477115e-18, -9.85508151674506e-19, -1.61235595962485e-18],
    [0, 0, 0, 0, 1.75775880167518e-18, 4.74982617364434e-19, -3.55429169456379e-19, 1.31185675272082e-18],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]])
# componentwise: 2 * abs_diff / (|Jacobian_1| + |Jacobian_2|)
rel_diff = np.array([
    [0, 0, 0, 0, 7.85175316554661e-08, -1.17485489544093e-07, 1.86115639593274e-08, 2.1131877425367e-07],
    [0, 0, 0, 0, -1.76422116783299e-08, 9.5190398570701e-08, 1.74871166493748e-07, 1.36484791784253e-07],
    [0, 0, 0, 0, 1.34933828136612e-07, 1.74871166493748e-07, -5.91304873754473e-08, -9.6741362219846e-08],
    [0, 0, 0, 0, 2.10931033342587e-07, 5.69979157415116e-08, -2.13257499690227e-08, 7.8711402111628e-08],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]])

# Masks: 0 ... tolerance met, 1 ... tolerance exceeded
abs_diff_mask = np.array([
    [0, 0, 0, 0, 1, 1, 0, 1],
    [0, 0, 0, 0, 0, 1, 1, 1],
    [0, 0, 0, 0, 1, 1, 0, 1],
    [0, 0, 0, 0, 1, 0, 0, 1],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]])
rel_diff_mask = np.array([
    [0, 0, 0, 0, 0, 1, 0, 1],
    [0, 0, 0, 0, 0, 0, 1, 1],
    [0, 0, 0, 0, 1, 1, 0, 0],
    [0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]])

M_1 = np.array([
    [1.08506944444444e-12, 5.42534722222222e-13, 2.71267361111111e-13, 5.42534722222222e-13, 0, 0, 0, 0],
    [5.42534722222222e-13, 1.08506944444444e-12, 5.42534722222222e-13, 2.71267361111111e-13, 0, 0, 0, 0],
    [2.71267361111111e-13, 5.42534722222222e-13, 1.08506944444444e-12, 5.42534722222222e-13, 0, 0, 0, 0],
    [5.42534722222222e-13, 2.71267361111111e-13, 5.42534722222222e-13, 1.08506944444444e-12, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]])
M_2 = np.array([
    [1.08506944444444e-12, 5.42534722222222e-13, 2.71267361111111e-13, 5.42534722222222e-13, 0, 0, 0, 0],
    [5.42534722222222e-13, 1.08506944444444e-12, 5.42534722222222e-13, 2.71267361111111e-13, 0, 0, 0, 0],
    [2.71267361111111e-13, 5.42534722222222e-13, 1.08506944444444e-12, 5.42534722222222e-13, 0, 0, 0, 0],
    [5.42534722222222e-13, 2.71267361111111e-13, 5.42534722222222e-13, 1.08506944444444e-12, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]])
K_1 = np.array([
    [1.83133333766667e-10, -3.32833334416667e-11, -9.15666668833333e-11, -5.82833334416667e-11, 0, 0, 0, 0],
    [-6.6616666775e-11, 2.164666671e-10, -4.1616666775e-11, -1.0823333355e-10, 0, 0, 0, 0],
    [-1.0823333355e-10, -4.1616666775e-11, 2.164666671e-10, -6.6616666775e-11, 0, 0, 0, 0],
    [-5.82833334416667e-11, -9.15666668833333e-11, -3.32833334416667e-11, 1.83133333766667e-10, 0, 0, 0, 0],
    [0, 0, 0, 0, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06],
    [0, 0, 0, 0, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06]])
K_2 = np.array([
    [1.83133333766667e-10, -3.32833334416667e-11, -9.15666668833333e-11, -5.82833334416667e-11, 0, 0, 0, 0],
    [-6.6616666775e-11, 2.164666671e-10, -4.1616666775e-11, -1.0823333355e-10, 0, 0, 0, 0],
    [-1.0823333355e-10, -4.1616666775e-11, 2.164666671e-10, -6.6616666775e-11, 0, 0, 0, 0],
    [-5.82833334416667e-11, -9.15666668833333e-11, -3.32833334416667e-11, 1.83133333766667e-10, 0, 0, 0, 0],
    [0, 0, 0, 0, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06, -3.33333333333333e-06],
    [0, 0, 0, 0, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06, -1.66666666666667e-06],
    [0, 0, 0, 0, -1.66666666666667e-06, -3.33333333333333e-06, -1.66666666666667e-06, 6.66666666666667e-06]])
b_1 = np.array([0, 0, 0, 0, 0, 0, 0, 0])
b_2 = np.array([0, 0, 0, 0, 0, 0, 0, 0])

### counter: 1 (end)

A convenient visualization of the results is possible, e.g., with Spyder (or any other Python IDE):
spyder

@chleh chleh force-pushed the chleh:jac-test branch from 459a83f to afa5d7b Oct 10, 2018

@chleh chleh added the please review label Oct 10, 2018

@chleh

This comment has been minimized.

Copy link
Collaborator Author

chleh commented Oct 16, 2018

I just changed the checks for the M, K and b matrices as requested. I hope the changed code works.
Regarding multiple values for the tolerances: E.g., for HM on a Quad8 that would imply 4*1 + 8*2 = 20 d.o.f., i.e., the Jacobian has 20^2 = 400 coefficients. Specifying 400 tolerance values seems pretty impractical to me.

<jacobian_assembler>
<type>CompareJacobians</type>
<jacobian_assembler>
<type>CentralDifferences</type>

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

I'd put an analytical jacobian as an example. Norbert has one.

This comment has been minimized.

@Scinopode

Scinopode Oct 17, 2018

<jacobian_assembler>
   <type>Analytical</type>
</jacobian_assembler>

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

OK

@@ -41,6 +42,8 @@ class AbstractJacobianAssembler
std::vector<double>& /*local_Jac_data*/,
LocalCoupledSolutions const& /*coupled_solutions*/)
{
// TODO make pure virtual.

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

Please add a short why, and maybe how (if not obvious).

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

Why ... because implementing this can be easily forgotten in derived classes. E.g., it's not implemented in the central differences assembler if I remember correctly.
How ... obviously: make it pure virtual and add the missing implementations.

Do you really want to have these comments in the code?

//! about exceeded tolerances and assembled local matrices.
std::ofstream _log_file;

//! Counter used for identifying blocks in the \c _log_file. Is be

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

Second sentence starts with "It is", probably...

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

Thanks.

namespace
{
//! Dumps a \c double value as a Python script snippet.
void dump_py(std::ofstream& fh, std::string const& var, double const val)

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

minor; an ostream would be more general...

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

fh << var << " = " << val << '\n';
}

//! Dumps an Eigen vector as a Python script snippet.

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

Dumps an arbitrary vector (of type Vec) as a Python script snippet.

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

dump_py_vec(fh, var, val);
}

template <typename Derived>

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

ah, ... and the Eigen vector comment goes here...

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator


if (fatal_error)
{
OGS_FATAL("%s", msg_fatal.c_str());

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

Not sure if the _log_file is flushed on std::abort (on throw it probably is correctly closed). I'd add a flush to the _log_file before every OGS_FATAL, just in case...

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const& config)
{
// TODO doc script corner case: Parameter could occur at different
// locations.

This comment has been minimized.

@endJunction

endJunction Oct 16, 2018

Member

note to myself. ?.

This comment has been minimized.

@chleh

chleh Oct 17, 2018

Author Collaborator

The parameter could be at
OpenGeoSys/processes/process/jacobian_assembler/type
or at
OpenGeoSys/processes/process/jacobian_assembler/CompareJacobians/jacobian_assembler/type
in the prj file. AFAIK this is not covered by the input parameter documentation scripts.

@endJunction
Copy link
Member

endJunction left a comment

Minor comments. Nice read!

@endJunction

This comment has been minimized.

Copy link
Member

endJunction commented Oct 16, 2018

@Scinopode If this works for you, give thumbs up, please.

@Scinopode

This comment has been minimized.

Copy link

Scinopode commented Oct 17, 2018

Works great, but output is still not complete (not flushed?)

@chleh chleh force-pushed the chleh:jac-test branch from 65a8664 to feb6e90 Oct 17, 2018

@wenqing
Copy link
Member

wenqing left a comment

Nice feature!

@endJunction endJunction merged commit c12d171 into ufz:master Oct 21, 2018

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/jenkins/pr-merge This commit looks good
Details
deploy/netlify Deploy preview ready!
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.