Peano 4
Loading...
Searching...
No Matches
Derivative.cpp
Go to the documentation of this file.
1#include "Derivative.h"
2#include "peano4/utils/Loop.h"
3
4#include <iostream>
5
6
8 int numberOfDoFsPerAxisInPatch,
9 const double* __restrict__ source,
10 int sourceIndex,
11 int sourceUnknowns,
12 int sourceHalo,
13 double* __restrict__ dest,
14 const tarch::la::Vector<Dimensions,int>& destIndex,
15 int destUnknowns,
16 int destHalo,
18) {
19 dfor(k,numberOfDoFsPerAxisInPatch) {
20 for (int d=0; d<Dimensions; d++) {
24
25 if (k(d)<=numberOfDoFsPerAxisInPatch/2) {
26 sourceVolumeRight(d)++;
27 }
28 else {
29 sourceVolumeLeft(d)--;
30 }
31
32 int sourceVolumeLeftLinearised = peano4::utils::dLinearised(sourceVolumeLeft, numberOfDoFsPerAxisInPatch+2*sourceHalo);
33 int sourceVolumeRightLinearised = peano4::utils::dLinearised(sourceVolumeRight, numberOfDoFsPerAxisInPatch+2*sourceHalo);
34 int destVolumeLinearised = peano4::utils::dLinearised(destVolume, numberOfDoFsPerAxisInPatch+2*destHalo);
35
36 double grad = ( source[sourceVolumeRightLinearised*sourceUnknowns+sourceIndex] - source[sourceVolumeLeftLinearised*sourceUnknowns+sourceIndex] ) / volumeH(d) / numberOfDoFsPerAxisInPatch;
37
38 dest[ destVolumeLinearised*destUnknowns+destIndex(d) ] = grad;
39 }
40 }
41}
42
43
45 int numberOfDoFsPerAxisInPatch,
46 const double* __restrict__ source,
47 int sourceIndex,
48 int sourceUnknowns,
49 int sourceHalo,
50 double* __restrict__ dest,
51 const tarch::la::Vector<Dimensions,int>& destIndex,
52 int destUnknowns,
53 int destHalo,
55) {
56 double result = 0.0;
57 dfor(k,numberOfDoFsPerAxisInPatch) {
58 for (int d=0; d<Dimensions; d++) {
62
63 if (k(d)<=numberOfDoFsPerAxisInPatch/2) {
64 sourceVolumeRight(d)++;
65 }
66 else {
67 sourceVolumeLeft(d)--;
68 }
69
70 int sourceVolumeLeftLinearised = peano4::utils::dLinearised(sourceVolumeLeft, numberOfDoFsPerAxisInPatch+2*sourceHalo);
71 int sourceVolumeRightLinearised = peano4::utils::dLinearised(sourceVolumeRight, numberOfDoFsPerAxisInPatch+2*sourceHalo);
72 int destVolumeLinearised = peano4::utils::dLinearised(destVolume, numberOfDoFsPerAxisInPatch+2*destHalo);
73
74 double grad = ( source[sourceVolumeRightLinearised*sourceUnknowns+sourceIndex] - source[sourceVolumeLeftLinearised*sourceUnknowns+sourceIndex] ) / volumeH(d) / numberOfDoFsPerAxisInPatch;
75
76 double delta = std::abs( dest[ destVolumeLinearised*destUnknowns+destIndex(d) ] - grad );
77 dest[ destVolumeLinearised*destUnknowns+destIndex(d) ] = grad;
78 result = std::max(result,delta);
79 }
80 }
81 return result;
82}
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define dfor(counter, max)
d-dimensional Loop
Definition Loop.h:308
CPUGPUMethod int dLinearised(const tarch::la::Vector< Dimensions, int > &counter, int max)
Map d-dimensional vector onto integer index.
Definition Loop.cpp:106
void computeGradient(int numberOfDoFsPerAxisInPatch, const double *__restrict__ source, int sourceIndex, int sourceUnknowns, int sourceHalo, double *__restrict__ dest, const tarch::la::Vector< Dimensions, int > &destIndex, int destUnknowns, int destHalo, const tarch::la::Vector< Dimensions, double > &volumeH)
This routine assumes that we have two patches of the same numberOfDoFsPerAxisInPatch.
Definition Derivative.cpp:7
double computeGradientAndReturnMaxDifference(int numberOfDoFsPerAxisInPatch, const double *__restrict__ source, int sourceIndex, int sourceUnknowns, int sourceHalo, double *__restrict__ dest, const tarch::la::Vector< Dimensions, int > &destIndex, int destUnknowns, int destHalo, const tarch::la::Vector< Dimensions, double > &volumeH)
Simple vector class.
Definition Vector.h:134