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.

Sign up for GitHub

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

Jump to bottom

Jacobian tester for local assemblers. #2238

Merged
merged 12 commits into from
Oct 21, 2018
Merged

Jacobian tester for local assemblers. #2238

merged 12 commits into from
Oct 21, 2018

Conversation

chleh
Copy link
Collaborator

@chleh 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
Copy link
Collaborator Author

chleh commented Oct 16, 2018 •
edited
Loading

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.

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
.../ProjectFile/prj/processes/process/jacobian_assembler/CompareJacobians/c_CompareJacobians.md Outdated
<jacobian_assembler>
<type>CompareJacobians</type>
<jacobian_assembler>
<type>CentralDifferences</type>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/AbstractJacobianAssembler.h
@@ -41,6 +42,8 @@ class AbstractJacobianAssembler
std::vector<double>& /*local_Jac_data*/,
LocalCoupledSolutions const& /*coupled_solutions*/)
{
// TODO make pure virtual.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/CompareJacobiansJacobianAssembler.h Outdated
//! about exceeded tolerances and assembled local matrices.
std::ofstream _log_file;

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks.

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/CompareJacobiansJacobianAssembler.cpp Outdated
namespace
{
//! Dumps a \c double value as a Python script snippet.
void dump_py(std::ofstream& fh, std::string const& var, double const val)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

minor; an ostream would be more general...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✅

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/CompareJacobiansJacobianAssembler.cpp Outdated
fh << var << " = " << val << '\n';
}

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✅

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/CompareJacobiansJacobianAssembler.cpp
dump_py_vec(fh, var, val);
}

template <typename Derived>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✅

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/CompareJacobiansJacobianAssembler.cpp

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

✅

endJunction
endJunction reviewed Oct 16, 2018
View reviewed changes
ProcessLib/CompareJacobiansJacobianAssembler.cpp
createCompareJacobiansJacobianAssembler(BaseLib::ConfigTree const& config)
{
// TODO doc script corner case: Parameter could occur at different
// locations.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

note to myself. ?.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
endJunction approved these changes Oct 16, 2018
View reviewed changes
Copy link
Member

@endJunction endJunction left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor comments. Nice read!

@endJunction
Copy link
Member

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

@Scinopode
Copy link
Contributor

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

wenqing
wenqing approved these changes Oct 18, 2018
View reviewed changes
Copy link
Member

@wenqing wenqing left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice feature!

@endJunction endJunction merged commit c12d171 into ufz:master Oct 21, 2018
@ogsbot
Copy link
Member

ogsbot commented Jun 19, 2020

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@Scinopode Scinopode Scinopode left review comments

@wenqing wenqing wenqing approved these changes

@endJunction endJunction endJunction approved these changes

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
5 participants
@chleh @endJunction @Scinopode @ogsbot @wenqing