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

Volumetric source terms #2220

Merged
merged 11 commits into from Sep 27, 2018

Conversation

Projects
None yet
3 participants
@TomFischer
Copy link
Member

TomFischer commented Sep 25, 2018

First version of volumetric source terms. Source term function has to be defined on the entire domain.

@chleh

chleh approved these changes Sep 25, 2018

Copy link
Collaborator

chleh left a comment

Only minor things.

{
DBUG("Assemble NodalSourceTerm.");
DBUG("Integrate NodalSourceTerm.");

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

Hmm... A nodal source term is not really integrated...

This comment has been minimized.

@TomFischer

TomFischer Sep 25, 2018

Author Member

Renamed.

class SourceTerm
{
public:
SourceTerm(const NumLib::LocalToGlobalIndexMap& source_term_dof_table)

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

Minor: explicit.

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

some parts are lost, it seems... But when you add mesh and variable's id to it, it's not needed then...

This comment has been minimized.

@endJunction

endJunction Sep 26, 2018

Member

explicit still missing...


#include "VolumetricSourceTerm.h"

#include <cassert>

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

not needed I guess.

: volumetric_source_term(volumetric_source_term_)
{}

VolumetricSourceTermData(VolumetricSourceTermData&& other)

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

Maybe now we don't have to specify that explicitly anymore due to newer compilers.


struct VolumetricSourceTermData final
{
VolumetricSourceTermData(

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

Minor: explicit

TESTER vtkdiff
REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
DIFF_DATA
square_1x1_quad_1e2_volumetricsourceterm_analytical_solution.vtu square_1e2_volumetricsourceterm_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure 1e-14 1e-16

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

Is the solution a linear slope or a constant? OGS is extremely close to the analytical result.

This comment has been minimized.

@TomFischer

TomFischer Sep 25, 2018

Author Member

The analytical solution is h(x,y) = - 1/2 (x^2 + x) + 2.

@@ -0,0 +1 @@
Declares a source term to be associated to the entire domain.

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

Maybe Declares a source term that is defined on the entire domain.

@@ -0,0 +1,2 @@
The name of the parameter that defines value that should be used for the source
term.

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

the value.

Maybe a comment that nodal parameters are not supported would be appropriate. And maybe there should also be an addendum that the user should carefully check which physical quantity that parameter is. It's not necessarily the same as \dot{primary_variable}.

This comment has been minimized.

@TomFischer

TomFischer Sep 25, 2018

Author Member

I do not understand what is meant by 'nodal parameter'. Can you explain, please?

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

Suggestion to keep the comment as is, and bake another PR handling this case.


The solution of this problem is
$$
h(x,y) = - 1 \over 2 (x^2 + x) + 2.

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

the fraction seems to be set wrongly. Btw. in LaTex there is \frac{...}{...}.

[menu.benchmarks]
parent = "elliptic"

+++

This comment has been minimized.

@chleh

chleh Sep 25, 2018

Collaborator

The test is OK and nicely documented, but for an actual user it might be more instructive to have some physical quantities with some actual units.

This comment has been minimized.

@TomFischer

TomFischer Sep 25, 2018

Author Member

You are right. I'll ask @waltherm or @Thomas-TK for a more practical test.

#pragma once

#include <memory>
#include "ProcessLib/Process.h"

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

Use forward decls like in CreateBoundaryCondition.h ...


private:
NumLib::LocalToGlobalIndexMap const& _dof_table;
std::size_t const _bulk_mesh_id;
MeshLib::Mesh const& _st_mesh;
int const _variable_id;

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

I'd move both, the variable's id and the mesh to the base class. The mesh is the domain of definition.
The variable's id is also part of domain of definition.

This comment has been minimized.

@TomFischer

TomFischer Sep 25, 2018

Author Member

I also thought about this. The VolumetricSourceTerm is also derived from the class SourceTerm and it doesn't need the mesh and the variable id as attributes. In the VolumetricSourceTerm implementation the mesh and the variable id are only used in the constructor.

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

OK, got it. Keep it as is.

source_term_dof_table.getNumberOfVariableComponents(variable_id))
{
OGS_FATAL(
"Variable id or component id too high. Actual values: (%d, %d), "

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

Could be split in two errors...


#include "SourceTerm.h"
#include "VolumetricSourceTermFEM.h"
#include "VolumetricSourceTermData.h"

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

Some of the includes could be replaced with fwd decls.

VolumetricSourceTermData(VolumetricSourceTermData&&) = default;

//! Copies are forbidden.
VolumetricSourceTermData(VolumetricSourceTermData const&) = delete;

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

Question: If the base class has a deleted assignment op or ctor, is a default assignment op or ctor created implicitly in the derived class?

This comment has been minimized.

@TomFischer

TomFischer Sep 26, 2018

Author Member

cppreference:
A defaulted copy assignment operator for class T is defined as deleted if T has a virtual base class that cannot be copy-assigned (selects a deleted or inaccessible function).

I'm unsure if this is an answer to your question?

### Comparison of the analytical solution and the computed solution

{{< img src="../square_1e2_volumetricsourceterm_pcs_0_ts_1_t_1.000000_h_AnalyticalSolution_VolumetricSourceTerm.png" >}}
{{< img src="../square_1e2_volumetricsourceterm_pcs_0_ts_1_t_1.000000_Diff_h_AnalyticalSolution_VolumetricSourceTerm.png" >}}

This comment has been minimized.

@endJunction

endJunction Sep 25, 2018

Member

Description of both figures is missing. The analytical solution line is not visible. In the second figure h should be pressure, I guess.

This comment has been minimized.

@TomFischer

TomFischer Sep 25, 2018

Author Member

In the formulae above the primary variable was 'h'. I renamed h -> p to make it consistent.

@TomFischer TomFischer force-pushed the TomFischer:VolumetricSourceTerms branch from 63c3fab to 3e91419 Sep 26, 2018

#include <vector>

#include "SourceTerm.h"
#include "VolumetricSourceTerm.h"

This comment has been minimized.

@endJunction

endJunction Sep 26, 2018

Member

self include? Or am I missing the difference?

virtual void integrate(const double t, GlobalVector& b) const = 0;

virtual ~SourceTerm() = default;
protected:

This comment has been minimized.

@endJunction

endJunction Sep 26, 2018

Member

Why protected? Could be private as far as I see. It can only be set in the ctor...

This comment has been minimized.

@TomFischer

TomFischer Sep 27, 2018

Author Member

Protected because the attribute should be used in the derived classes.


private:
Parameter<double> const& _volumetric_source_term;
std::vector<std::unique_ptr<VolumetricSourceTermLocalAssemblerInterface>>

This comment has been minimized.

@endJunction

endJunction Sep 26, 2018

Member

I'd think a fwd-decl to the local asm interface would be sufficient...

This comment has been minimized.

@TomFischer

TomFischer Sep 27, 2018

Author Member

The fwd-decl leads to:

error: invalid application of 'sizeof' to an incomplete type 'ProcessLib::VolumetricSourceTermLocalAssemblerInterface'

template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
class VolumetricSourceTermLocalAssembler

This comment has been minimized.

@endJunction

endJunction Sep 26, 2018

Member

final

auto const st_val = _volumetric_source_term(t, pos)[0];

_local_rhs.noalias() +=
sm.N * st_val * sm.detJ * sm.integralMeasure * wp.getWeight();

This comment has been minimized.

@endJunction

endJunction Sep 26, 2018

Member

It would be much better to store only the N vector in the ip data (not all of the shape matrices) and precompute the integration_weight = detJ * intMeasure * w_ip.

This comment has been minimized.

@TomFischer

TomFischer Sep 27, 2018

Author Member

Introduced class IntegrationPointData, see commit 3ba3927.

@endJunction
Copy link
Member

endJunction left a comment

Looks good in general but still few things either are not yet uploaded, or forgotten.

@TomFischer TomFischer force-pushed the TomFischer:VolumetricSourceTerms branch from 3e91419 to 3ba3927 Sep 27, 2018

@TomFischer TomFischer force-pushed the TomFischer:VolumetricSourceTerms branch from 3ba3927 to 79c3e71 Sep 27, 2018

@endJunction endJunction merged commit 6ec6371 into ufz:master Sep 27, 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

@TomFischer TomFischer deleted the TomFischer:VolumetricSourceTerms branch Sep 27, 2018

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.