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

Update couple your code #315

Merged
merged 3 commits into from
Dec 22, 2023
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -15,30 +15,36 @@ For volume coupling in 2D, mesh connectivity boils down to defining triangles an

All kind of connectivity can be built up directly from vertices. Triangles and quads also allow us to define them using edge IDs.

<!-- TODO: What about setMeshEdges, setMeshTriangle, setMeshQuad, setMeshTetrahedron? -->
```cpp
int setMeshEdge (precice::string_view meshName, int firstVertexID, int secondVertexID);
void setMeshTriangle (precice::string_view meshName, int firstVertexID, int secondVertexID, int thirdVertexID);
void setMeshQuad(precice::string_view meshName, int firstVertexID, int secondVertexID, int thirdVertexID, int fourthVertexID);
void setMeshTetrahedron(precice::string_view meshName, int firstVertexID, int secondVertexID, int thirdVertexID, int fourthVertexID);
void setMeshEdge(precice::string_view meshName, VertexID firstVertexID, VertexID secondVertexID);
void setMeshTriangle(precice::string_view meshName, VertexID firstVertexID, VertexID secondVertexID, VertexID thirdVertexID);
void setMeshQuad(precice::string_view meshName, VertexID firstVertexID, VertexID secondVertexID, VertexID thirdVertexID, VertexID fourthVertexID);
void setMeshTetrahedron(precice::string_view meshName, VertexID firstVertexID, VertexID secondVertexID, VertexID thirdVertexID, VertexID fourthVertexID);
```

* `setMeshEdge` defines a mesh edge between two vertices and returns an edge ID.
* `setMeshTriangle` defines a mesh triangle by three edges.
* `setMeshQuad` defines a mesh quad by four edges.
* `setMeshTetrahredron` defines a mesh tetrahedron by four vertices.
There are also bulk versions of these methods, which can be easier to handle in some cases:

```cpp
void setMeshEdges(precice::string_view meshName, precice::span<const VertexID> vertices);
void setMeshTriangles(precice::string_view meshName, precice::span<const VertexID> vertices);
void setMeshQuads(precice::string_view meshName, precice::span<const VertexID> vertices);
void setMeshTetrahedra(precice::string_view meshName, precice::span<const VertexID> vertices);
```

If you do not configure any features in the preCICE configuration that require mesh connectivity, all these API functions are [no-ops](https://en.wikipedia.org/wiki/NOP_(code)). Thus, don't worry about performance. If you need a significant workload to already create this connectivity information in your adapter in the first place, you can also explicitly ask preCICE whether it is required:

<!-- TODO: needs update -->
```cpp
bool isMeshConnectivityRequired(int meshID);
bool requiresMeshConnectivityFor(precice::string_view meshName);
```

{% warning %}
The API function `isMeshConnectivityRequired` is only supported since v2.3.
{% endwarning %}

{% warning %}
The bulk API functions are only supported since v3.
fsimonis marked this conversation as resolved.
Show resolved Hide resolved
{% endwarning %}

Maybe interesting to know: preCICE actually does internally not compute with quads, but creates two triangles. [Read more](https://precice.discourse.group/t/highlights-of-the-new-precice-release-v2-1/274#2-1-using-quads-for-projection).

{% warning %}
Expand All @@ -48,52 +54,45 @@ Quads are only supported since v2.1. For older version, the methods only exist a
The following code shows how mesh connectivity can be defined in our example. For sake of simplification, let's only define one triangle and let's assume that it consists of the first three vertices.

```cpp

[...]

std::vector<int> vertexIDs(vertexSize);
/* ... */

// We define the unit square in 2D
std::vector<double> coords{
0, 0, // A
0, 1, // B
1, 0, // C
1, 1 // D
};
std::vector<VertexID> vertexIDs(4);
IshaanDesai marked this conversation as resolved.
Show resolved Hide resolved
precice.setMeshVertices(meshName, coords, vertexIDs);

int edgeIDs[3];
edgeIDs[0] = precice.setMeshEdge(meshName, vertexIDs[0], vertexIDs[1]);
edgeIDs[1] = precice.setMeshEdge(meshName, vertexIDs[1], vertexIDs[2]);
edgeIDs[2] = precice.setMeshEdge(meshName, vertexIDs[2], vertexIDs[0]);

if(dim==3)
precice.setMeshTriangle(meshName, edgeIDs[0], edgeIDs[1], edgeIDs[2]);

[...]
if (precice.requiresMeshConnectivityFor(meshName)) {

// defines triangles ABC and BCD separately
precice.setMeshTriangle(meshName, vertexIDs[0], vertexIDs[1], vertexIDs[2]);
precice.setMeshTriangle(meshName, vertexIDs[1], vertexIDs[2], vertexIDs[3]);

// defines triangles ABC and BCD in one go
std::vector<VertexID> triangles{
vertexIDs[0], vertexIDs[1], vertexIDs[2],
vertexIDs[1], vertexIDs[2], vertexIDs[3]
};
precice.setMeshTriangles(meshName, triangles)
}

/* ... */
```

<!-- TODO: Section below should be the default on branch precice-v3 -->
## Changes in v3

Version 3 overhauls the definition of meshes.

Connectivity now consists of explicitly defined elements (elements created via calls to the API) and implicitly defined elements (elements additionally created by preCICE).
As an example, explicitly defining a triangle ABC via the API guarantees the existence of the implicit edges AB, AC, and BC.

Furthermore, all connectivity is defined only via vertex IDs. There are no more edge IDs to worry about.
The order of vertices also does not matter. Triangles BAC and CAB are considered duplicates and preCICE removes one of them during the deduplication step.
## Mesh pre-processing

The API for defining individual connectivity elements looks as follows:
preCICE pre-processes all provided meshes during initialization, removing duplicates and adding missing hierarchical elements.

```cpp
void setMeshEdge(precice::string_view meshName, int firstVertexID, int secondVertexID);
void setMeshTriangle(precice::string_view meshName, int firstVertexID, int secondVertexID, int thirdVertexID);
void setMeshQuad(precice::string_view meshName, int firstVertexID, int secondVertexID, int thirdVertexID, int fourthVertexID);
void setMeshTetrahredron(precice::string_view meshName, int firstVertexID, int secondVertexID, int thirdVertexID, int fourthVertexID);
```

Each of the above functions is accompanied by a bulk version, which allows to set multiple elements in a single call.
Some projection-based features require all hierarchical elements to be present.
Meaning, a triangle ABC requires edges AB, BC and AC to exist.
Manually defining such elements is not a pleasant experience, especially when dealing with tetrahedra whilst avoiding duplicates.

```cpp
void setMeshEdges(precice::string_view meshName, ::precice::span<const VertexID> vertices);
void setMeshTriangles(precice::string_view meshName, ::precice::span<const VertexID> vertices);
void setMeshQuads(precice::string_view meshName, ::precice::span<const VertexID> vertices);
void setMeshTetrahedra(precice::string_view meshName, ::precice::span<const VertexID> vertices);
```
This is why preCICE steps in and handles this internally.
In practise, you only need to define the connectivity your solver exposes.

## Putting it all together

Expand Down
49 changes: 28 additions & 21 deletions pages/docs/couple-your-code/couple-your-code-direct-access.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,41 +14,48 @@ These API functions are work in progress, experimental, and are not yet released
This concept is required if you want to access received meshes directly. It might be relevant in case you don't want to use the mapping schemes in preCICE, but rather want to use your own solver for data mapping. As opposed to the usual preCICE mapping, only a single mesh (from the other participant) is now involved in this situation since an 'own' mesh defined by the participant itself is not required any more. In order to re-partition the received mesh, the participant needs to define the mesh region it wants read data from and write data to. The complete concept on the receiving participant looks as follows:

```cpp
// Allocate a bounding-box vector containing lower and upper bounds per
// space dimension
std::vector<double> boundingBox(dim * 2);

// fill the allocated 'boundingBox' according to the interested region
// with the desired bounds...
// Get relevant IDs. Note that "ReceivedMeshname" is not a name of a
// Note that "ReceivedMeshname" is not a name of a
// provided mesh, but a mesh defined by another participant. Accessing
// a received mesh directly is disabled in a usual preCICE configuration.
const int otherMeshID = precice.getMeshID("ReceivedMeshName");
const int writeDataID = precice.getDataID("WriteDataName", otherMeshID);
const std::string otherMesh = "ReceivedMeshName";

// Get the spacial dimensionality of the mesh
const int dims = precice.getMeshDimensions(otherMesh);

// Allocate and fill the 'boundingBox' according to the interested region
// with the desired bounds, in our example we use the unit cube.
// Assuming dim == 3, means that the bounding box has dim * 2 == 6 elements.
std::vector<double> boundingBox {
0, 0, 0, // minimum corner
1, 1, 1, // maximum corner
fsimonis marked this conversation as resolved.
Show resolved Hide resolved
};

// Define region of interest, where we want to obtain the direct access.
// See also the API documentation of this function for further notes.
precice.setMeshAccessRegion(otherMeshID, boundingBox.data());
precice.setMeshAccessRegion(otherMesh, boundingBox);

// initialize preCICE as usual
double dt = precice.initialize();
precice.initialize();

// Get the size of the received mesh partition, which lies within the
// defined bounding (provided by the coupling participant)
const int otherMeshSize = precice.getMeshVertexSize(otherMeshID);
const int otherMeshSize = precice.getMeshVertexSize(otherMesh);

// Now finally get the data. First allocate memory for the IDs and the
// vertices
std::vector<double> otherSolverVertices(otherMeshSize * dim);
std::vector<int> ids(otherMeshSize);
// Now finally get information about the mesh vertices.
// First allocate memory for the IDs and coordinates
std::vector<double> otherCoordinates(otherMeshSize * dim);
fsimonis marked this conversation as resolved.
Show resolved Hide resolved
std::vector<VertexID> otherVertexIDs(otherMeshSize);
// ... and afterwards ask preCICE to fill the vectors
precice.getMeshVerticesAndIDs(otherMeshID,
otherMeshSize,
ids.data(),
otherSolverVertices.data());
precice.getMeshVertexIDsAndCoordinates(otherMesh,
otherVertexIDs,
otherCoordinates);

// continue with time loop and write data directly using writeDataID and
// continue with time loop and write data directly to the mesh using
// the received ids, which correspond to the vertices
const int dataDims = precice.getDataDimensions(otherMesh, "OtherData");
std::vector<double> data(dataDims * otherMeshSize);
precice.writeData(otherMesh, "OtherData", otherVertexIDs, data);

```

## Concept and API
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,12 @@ keywords: api, adapter, parallelization, mpi, initialization
summary: "preCICE uses MPI for communication between different participants (and also for communication between ranks of the same participant). So are there any problems if the solver that you intend to couple also already uses MPI (e.g. for parallelization)? Who should initialize MPI? Who should finalize MPI? This is what we discuss here."
---

It is not complicated. There are just three rules that preCICE follows:
It is not complicated. There are three rules that preCICE follows:

* preCICE only initializes MPI if it is not yet initialized (by e.g. the solver you want to couple).
* preCICE finalizes MPI if and only if it was also initialized by preCICE.
* preCICE only initializes MPI if it needs MPI.
* preCICE uses `MPI_COMM_WORLD` if no custom communicator is provided.

So what does this mean for your adapter code:

Expand All @@ -26,8 +27,7 @@ MPI_Comm_size(MPI_COMM_WORLD, &world_size);

[...] // maybe more initialization

precice::Participant precice("SolverName", world_rank, world_size);
precice.configure("precice-config.xml");
precice::Participant precice("SolverName", "precice-config.xml", world_rank, world_size);

[...] // declare meshes vertices etc.

Expand All @@ -42,3 +42,17 @@ precice.finalize();
MPI_Finalize();

```

If you need to provide a custom communicator to preCICE, you can do so as follows:

```cpp

[...] // Initialize MPI and start your solver

MPI_Comm my_comm = /* ... */;

precice::Participant precice("SolverName", "precice-config.xml", world_rank, world_size, &my_comm);

[...]

```
8 changes: 4 additions & 4 deletions pages/docs/couple-your-code/couple-your-code-gradient-data.md
Original file line number Diff line number Diff line change
Expand Up @@ -36,14 +36,14 @@ Let's add gradient data to our example code:
precice::Participant precice("FluidSolver", "precice-config.xml", rank, size); // constructor

int dim = precice.getMeshDimensions("FluidMesh");
[...]
/* ... */
precice.setMeshVertices("FluidMesh", vertexSize, coords, vertexIDs);

std::vector<double> stress(vertexSize * dim);

// create gradient data
std::vector<double> stressGradient(vertexSize * dim * dim)
[...]
/* ... */
precice.initialize();

while (not simulationDone()){ // time loop
Expand All @@ -58,14 +58,14 @@ while (not simulationDone()){ // time loop
precice.writeData("FluidMesh", "Stress", vertexIDs, stress);

// write gradient data
if (isGradientDataRequired(dataID)){
if (requiresGradientDataFor("FluidMesh")){
computeStressGradient(stressGradient)
precice.writeGradientData("FluidMesh", "Stress", vertexIDs, stressGradient);
}

precice.advance(dt);
}
[...]
/* ... */
```

{% experimental %}
Expand Down
28 changes: 10 additions & 18 deletions pages/docs/couple-your-code/couple-your-code-implicit-coupling.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,17 @@ keywords: api, adapter, coupling schemes, checkpoint, fixed-point
summary: "In previous steps, we only considered explicit coupling. We now move onto implicit coupling, so sub-iterating each time step multiple times until a convergence threshold is reached. This stabilzes strongly-coupled problems."
---

The main ingredient needed for implicit coupling is move backwards in time. For that, we need a [flux capacitor](https://www.youtube.com/watch?v=VcZe8_RZO8c). Just kidding :wink:. What we really need is that your solver can write and read iteration checkpoints. An iteration checkpoint should contain all the information necessary to reload a previous state of your solver. What exactly is needed depends solely on your solver. preCICE tells you when you need to write and read checkpoints. To this end, preCICE uses the following action interface:
The main ingredient needed for implicit coupling is move backwards in time. For that, we need a [flux capacitor](https://www.youtube.com/watch?v=VcZe8_RZO8c). Just kidding :wink:. What we really need is that your solver can write and read iteration checkpoints. An iteration checkpoint should contain all the information necessary to reload a previous state of your solver. What exactly is needed depends solely on your solver. preCICE tells you when you need to write and read checkpoints. To this end, preCICE uses the following interface:

```cpp
bool isActionRequired(const std::string& action)
void markActionFulfilled(const std::string& action)
const std::string& constants::actionReadIterationCheckpoint()
const std::string& constants::actionWriteIterationCheckpoint()
bool requiresWritingCheckpoint()
bool requiresReadingCheckpoint()
```

* `isActionRequired` inquires the necessity of a certain action. It takes a string argument to reference the action.
* `markActionFulfilled` tells preCICE that the action is fulfilled. This is a simple safeguard. If a certain action is required and you did not mark it as fulfilled preCICE will complain.
* The Methods in the `precice::constants` namespace return strings to reference specific actions. For implicit coupling, we need `actionReadIterationCheckpoint` and `actionWriteIterationCheckpoint`.
These functions perform double duty:

1. They inform the adapter that writing or reading a checkpoint is required by the solver.
2. They let preCICE know that your adapter is capable of implicit coupling. preCICE will show an error if you configure implicit coupling without calling these functions.

Let's extend our example code to also handle implicit coupling.

Expand All @@ -25,17 +24,13 @@ turnOnSolver(); //e.g. setup and partition mesh

precice::Participant precice("FluidSolver","precice-config.xml",rank,size); // constructor

const std::string& coric = precice::constants::actionReadIterationCheckpoint();
const std::string& cowic = precice::constants::actionWriteIterationCheckpoint();

int dim = precice.getMeshDimensions("FluidMesh");
int vertexSize; // number of vertices at wet surface
// determine vertexSize
std::vector<double> coords(vertexSize*dim); // coords of vertices at wet surface
// determine coordinates
std::vector<int> vertexIDs(vertexSize);
precice.setMeshVertices("FluidMesh", coords, vertexIDs);
delete[] coords;

std::vector<double> forces(vertexSize*dim);
std::vector<double> displacements(vertexSize*dim);
Expand All @@ -48,9 +43,8 @@ double dt; // actual time step size
```cpp
precice.initialize();
while (precice.isCouplingOngoing()){
if(precice.isActionRequired(cowic)){
if(precice.requiresWritingCheckpoint()){ // new time window
saveOldState(); // save checkpoint
precice.markActionFulfilled(cowic);
}
preciceDt = precice.getMaxTimeStepSize();
solverDt = beginTimeStep(); // e.g. compute adaptive dt
Expand All @@ -61,16 +55,14 @@ while (precice.isCouplingOngoing()){
computeForces(forces);
precice.writeData("FluidMesh", "Forces", vertexIDs, forces);
precice.advance(dt);
if(precice.isActionRequired(coric)){ // time step not converged
if(precice.requiresReadingCheckpoint()){ // iteration did not converge
reloadOldState(); // set variables back to checkpoint
precice.markActionFulfilled(coric);
}
else{ // time step converged
else{ // iteration converged
endTimeStep(); // e.g. update variables, increment time
}
}
precice.finalize(); // frees data structures and closes communication channels
delete[] vertexIDs, forces, displacements;
turnOffSolver();
```

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,6 @@ std::vector<double> coords(vertexSize*dim); // coords of vertices at wet surface
// determine coordinates
std::vector<int> vertexIDs(vertexSize);
precice.setMeshVertices("FluidMesh", coords, vertexIDs);
delete[] coords;

std::vector<double> forces(vertexSize*dim);
std::vector<double> displacements(vertexSize*dim);
Expand All @@ -73,7 +72,6 @@ while (not simulationDone()){ // time loop
endTimeStep(); // e.g. update variables, increment time
}
precice.finalize(); // frees data structures and closes communication channels
delete[] vertexIDs, forces, displacements;
turnOffSolver();
```

Expand Down
Loading