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

Balance for HT process #2132

Merged
merged 13 commits into from May 31, 2018

Conversation

Projects
None yet
3 participants
@TomFischer
Copy link
Member

TomFischer commented May 28, 2018

PR implements getFlux in the in the class HTProcess and the local assembler for the HT process. The new feature is checked by two new test cases.

@TomFischer TomFischer force-pushed the TomFischer:BalanceForHTProcess branch from 829473b to 100c15f May 28, 2018

@TomFischer TomFischer added the WIP 🏗 label May 28, 2018

@TomFischer TomFischer force-pushed the TomFischer:BalanceForHTProcess branch 2 times, most recently from 3ab0a3c to d64bad6 May 28, 2018

@TomFischer TomFischer removed the WIP 🏗 label May 28, 2018

fe.computeShapeFunctions(pnt_local_coords.getCoords(), shape_matrices,
GlobalDim, false);
std::vector<double> flux;
flux.resize(3);

This comment has been minimized.

@endJunction

endJunction May 29, 2018

Member

One could move it down to near the return statement.
I'd say, that returning a Eigen::Vector3d is safer, then the std::vector...

NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N, T_int_pt, p_int_pt);

// set the values into array vars
vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] =

This comment has been minimized.

@endJunction

endJunction May 29, 2018

Member

Would it make sense to pass the vars[T] directly to shapeFunctionInterpolate call?

NumLib::shapeFunctionInterpolate(local_x, shape_matrices.N,
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)],
    vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)]);
MaterialLib::Fluid::FluidPropertyType::Density, vars);
auto const b = this->_material_properties.specific_body_force;
// here it is assumed that the vector b is directed 'downwards'
q += K_over_mu * rho_w * b;

This comment has been minimized.

@endJunction

endJunction May 29, 2018

Member

I'm not getting the intentions of the comment. What is 'downwards'? Maybe it's copied from some other context...

DBUG("This is the thermal part of the staggered HTProcess.");
return;
}
if (_balance_mesh) // computing the balance is optional

This comment has been minimized.

@endJunction

endJunction May 29, 2018

Member
if (!_balance_mesh)
{
    return;
}
auto* const balance_pv = MeshLib::getOrCreateMeshProperty<double>(
*_balance_mesh, _balance_pv_name, MeshLib::MeshItemType::Cell, 1);
// initialise the PropertyVector pv with zero values
std::fill(balance_pv->begin(), balance_pv->end(), 0.0);

This comment has been minimized.

@endJunction

endJunction May 29, 2018

Member

Not sure here;
If the properties are set, then use NaN for initialization. If the properties are '+='-ed, then setting them 0 is of course necessary.

This comment has been minimized.

@TomFischer

TomFischer May 30, 2018

Author Member

The properties are '+='-ed.

@endJunction
Copy link
Member

endJunction left a comment

Smaller comments, mostly for discussion.
Didn't look into the ctests.

@TomFischer TomFischer force-pushed the TomFischer:BalanceForHTProcess branch from 59a9966 to 3901500 May 30, 2018

@chleh
Copy link
Collaborator

chleh left a comment

Two small comments and two bigger ones:

  • time dependence?
  • mass vs. heat flux?
@@ -105,7 +105,7 @@ class LocalAssemblerData : public GroundwaterFlowLocalAssemblerInterface
/// Computes the flux in the point \c p_local_coords that is given in local
/// coordinates using the values from \c local_x.
// TODO add time dependency
std::vector<double> getFlux(

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

Isn't the flux 2D in two dimentsions?

@@ -88,6 +88,73 @@ class HTFEM : public HTLocalAssemblerInterface
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}

/// Computes the flux in the point \c pnt_local_coords that is given in
/// local coordinates using the values from \c local_x.
// TODO add time dependency

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

Yep, that would be nice (for the sake of completeness). Could you please add it now?

auto const K =
this->_material_properties.porous_media_properties
.getIntrinsicPermeability(t, pos)
.getValue(t, pos, 0.0,

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

Is the 0.0 correct here?

This comment has been minimized.

@TomFischer

TomFischer May 30, 2018

Author Member

Yes, I think so.

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

OK, class Permeability says that this can have an arbitrary value.

}

Eigen::Vector3d flux;
flux.head(GlobalDim) = q;

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

Maybe here the templated version head>GlobalDim>() would be a tiny bit better.
And I see: That's why flux is always 3D.

This comment has been minimized.

@endJunction

endJunction May 30, 2018

Member

Regarding the flux being 3D always: I thought it would be easier in general to be it of constant size, like the node coordinates.

flux.head<GlobalDim>() would be better, that's true!

/// Computes the flux in the point \c pnt_local_coords that is given in
/// local coordinates using the values from \c local_x.
// TODO add time dependency
Eigen::Vector3d getFlux(

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

For the HT process the problem is that there are potentially two different interesting fluxes:
The mass/volume flux and the heat flux.
Maybe there should be different names for them?

This comment has been minimized.

@endJunction

endJunction May 30, 2018

Member

Or one could pass the variable_id to getFlux?

This comment has been minimized.

@TomFischer

TomFischer May 30, 2018

Author Member

You are right, there are potentially two different interesting fluxes. Maybe other processes have also interesting quantities fluxes. For this reason we would need rather general names. I opened an issue for discussion and for the time being I would keep the implementation of this PR. Is this okay for you?

This comment has been minimized.

@chleh

chleh May 30, 2018

Collaborator

OK.

@TomFischer TomFischer force-pushed the TomFischer:BalanceForHTProcess branch from 3901500 to a823565 May 30, 2018

@TomFischer

This comment has been minimized.

Copy link
Member Author

TomFischer commented May 30, 2018

Added the time to all getFlux and the balance functions, see commit a823565. For the second bigger comment I added the issue #2134 for discussion.

@chleh

chleh approved these changes May 30, 2018

Copy link
Collaborator

chleh left a comment

👍 thanks for the changes.

@TomFischer TomFischer force-pushed the TomFischer:BalanceForHTProcess branch from a823565 to dd8c5ca May 31, 2018

@TomFischer TomFischer merged commit f2c767a into ufz:master May 31, 2018

2 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/jenkins/pr-merge This commit looks good
Details

@TomFischer TomFischer deleted the TomFischer:BalanceForHTProcess branch May 31, 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.