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

3-Tensor Operations and Constructors #415

Merged
merged 17 commits into from
Oct 14, 2020
Merged

3-Tensor Operations and Constructors #415

merged 17 commits into from
Oct 14, 2020

Conversation

zjwegert
Copy link
Contributor

@zjwegert zjwegert commented Oct 7, 2020

Added functionality to compute

a_ij = b_kij*c_k
a_i = b_ijk*c_jk

Added functionality to compute
```
a_ij = b_kij*c_k
a_i = b_ijk*c_jk
```
@codecov-io
Copy link

codecov-io commented Oct 7, 2020

Codecov Report

Merging #415 into master will increase coverage by 0.05%.
The diff coverage is 93.50%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #415      +/-   ##
==========================================
+ Coverage   87.26%   87.32%   +0.05%     
==========================================
  Files         158      158              
  Lines       11094    11170      +76     
==========================================
+ Hits         9681     9754      +73     
- Misses       1413     1416       +3     
Impacted Files Coverage Δ
src/TensorValues/TensorValues.jl 100.00% <ø> (ø)
src/TensorValues/ThirdOrderTensorValueTypes.jl 52.77% <16.66%> (-7.23%) ⬇️
src/TensorValues/Operations.jl 87.24% <100.00%> (+4.13%) ⬆️
src/TensorValues/VectorValueTypes.jl 79.54% <0.00%> (-2.28%) ⬇️
src/ReferenceFEs/CLagrangianRefFEs.jl 94.52% <0.00%> (-0.31%) ⬇️
src/TensorValues/Indexing.jl 84.61% <0.00%> (+9.61%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4c94645...45ab2b3. Read the comment docs.

Copy link
Member

@fverdugo fverdugo left a comment

Choose a reason for hiding this comment

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

@Omega-xyZac thanks for your PR.

I have added some coments to be addressed (see inline).

In addition:

src/TensorValues/Operations.jl Outdated Show resolved Hide resolved
@zjwegert
Copy link
Contributor Author

zjwegert commented Oct 8, 2020

Bug reported in #414 is due to not using @law for part of the weak form. Unrelated to new functionality.

Omega-xyZac added 2 commits October 9, 2020 09:12
Forgot to switch inner and outer for loop
- This new functionality calculates `a_ijpm = b_ijkl*c_jkpm` and also provides testing.
- I have re-enabled use of `:` operator since this is the most suitable.
@zjwegert
Copy link
Contributor Author

zjwegert commented Oct 9, 2020

I guess we need a new symbol instead of : @fverdugo.

@fverdugo
Copy link
Member

fverdugo commented Oct 9, 2020

I guess we need a new symbol instead of : @fverdugo.

Sure @Omega-xyZac ! Perhaps we can use ⊡ (\boxdot)

In addition, try to respect our convention in which the indices that are contracted are the ones that are in "contact". That is:

a = b  c # a_ijkl = b_ijab*c_abkl

This convention is important since products with tensors are not always commutative.

@zjwegert
Copy link
Contributor Author

zjwegert commented Oct 9, 2020

I guess we need a new symbol instead of : @fverdugo.

Sure @Omega-xyZac ! Perhaps we can use ⊡ (\boxdot)

In addition, try to respect our convention in which the indices that are contracted are the ones that are in "contact". That is:

a = b  c # a_ijkl = b_ijab*c_abkl

This convention is important since products with tensors are not always commutative.

Ah yes sorry! That was actually a typo in the comment. The new \boxdot functionality adheres to that.

Omega-xyZac added 2 commits October 9, 2020 21:18
- Also changed comment to correctly reflect computation
@santiagobadia
Copy link
Member

Why do we introduce a new symbol for the double contraction of tensors? There is a standard symbol for this and it is : . I would keep : unless there is a very compelling reason.

@zjwegert
Copy link
Contributor Author

zjwegert commented Oct 9, 2020

Why do we introduce a new symbol for the double contraction of tensors? There is a standard symbol for this and it is : . I would keep : unless there is a very compelling reason.

I agree but there seems to be issues with :? Perhaps @fverdugo can explain this more.

@fverdugo
Copy link
Member

fverdugo commented Oct 9, 2020

@santiagobadia the symbol : has a very concrete meaning in julia (range). I find very confusing to give it another meaning.

@santiagobadia
Copy link
Member

There is not going to be any clash between these two cases. The problem is that the precedence of : vs +,- is different between these two cases. Probably, this cannot be fixed. That is a compelling reason. So 👍 with the new symbol.

@zjwegert
Copy link
Contributor Author

zjwegert commented Oct 9, 2020

Hi all,
I have some other contractions that I would like to add. One of these will break the convention a_ijmp = b_kij*c_kmp. I don't see any nice way around this. I can just leave this one out and use it in my personal code.

The others are

# a_ilm = b_ijk*c_jklm
# a_ilm = b_ij*c_jlm
# a_il = b_ijk*c_jkl
# a_il = b_ij*c_jl

@santiagobadia
Copy link
Member

The last operations are single and double contractions. You can add them and use the agreed symbols. I would not include the one that breaks the rules in Gridap. I don't think this is standard enough. You can define them in your part.

@santiagobadia
Copy link
Member

We could probably consider also this notation...

¹,²,³,⁴,... # n-contractions
¹ == 
² == 

I guess we will need 4-contractions of 4-tensors, which are usually expressed as :: or 3-contractions, etc.

What do you think @fverdugo @Omega-xyZac ?

@zjwegert
Copy link
Contributor Author

This is probably the best solution to the problem of higher order contractions.

@@ -220,8 +246,40 @@ function inner(a::SymFourthOrderTensorValue{D},b::MultiValue{Tuple{D,D}}) where
inner(a,symmetric_part(b))
end

# a_i = b_ijk*c_jk
@generated function inner(a::A, b::B) where {A<:MultiValue{Tuple{D1,D2,D3}},B<:MultiValue{Tuple{D2,D3}}} where {D1,D2,D3}
Copy link
Member

Choose a reason for hiding this comment

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

This is not an inner product. This is a double contraction, isn't it?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes. I was trying to follow the convention that has been used in the rest of the file. I'll update these with the above notation

###############################################################

# a_ijpm = b_ijkl*c_klpm
@generated function ⊡(a::A, b::B) where {A<:SymFourthOrderTensorValue{D},B<:SymFourthOrderTensorValue{D}} where D
Copy link
Member

Choose a reason for hiding this comment

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

I would give a name to this method, e.g., double_contraction and then define the symbol, as we do, e.g., for dot.

Omega-xyZac added 3 commits October 10, 2020 16:07
Added
```
# a_ilm = b_ijk*c_jklm
# a_ilm = b_ij*c_jlm
# a_il = b_ijk*c_jkl
# a_il = b_ij*c_jl
```
Implemented notation per S. Badia:
```
⋅¹,⋅²,⋅³,⋅⁴,... # n-contractions
⋅¹ == ⋅
⋅² == ⊡
```
Remove ⊡ from op
@zjwegert
Copy link
Contributor Author

Noticed this in the build doc

$ julia --check-bounds=yes --color=yes -e "if VERSION < v\"0.7.0-DEV.5183\"; Pkg.test(\"${JL_PKG}\", coverage=true); else using Pkg; Pkg.test(coverage=true); end"

   Testing Gridap

 Resolving package versions...

WARNING: Method definition dot(A<:(Gridap.TensorValues.MultiValue{Tuple{D1, D2}, T, N, L} where L where N where T), B<:(Gridap.TensorValues.MultiValue{Tuple{D2, D3}, T, N, L} where L where N where T)) where {D1, D2, D3, A<:(Gridap.TensorValues.MultiValue{Tuple{D1, D2}, T, N, L} where L where N where T), B<:(Gridap.TensorValues.MultiValue{Tuple{D2, D3}, T, N, L} where L where N where T)} in module TensorValues at /home/travis/build/gridap/Gridap.jl/src/TensorValues/Operations.jl:158 overwritten at /home/travis/build/gridap/Gridap.jl/src/TensorValues/Operations.jl:212.

  ** incremental compilation may be fatally broken for this module **

WARNING: Method definition dot(Gridap.Helpers.GridapType, Gridap.Helpers.GridapType) in module Helpers at /home/travis/build/gridap/Gridap.jl/src/Helpers/GridapTypes.jl:41 overwritten in module TensorValues at /home/travis/build/gridap/Gridap.jl/src/TensorValues/Operations.jl:600.

  ** incremental compilation may be fatally broken for this module **

I think this is due to lines 339,340, and 598-606. Is this the correct way to be doing this? Perhaps ⋅¹ should not be in like 598.

Omega-xyZac and others added 4 commits October 10, 2020 16:38
Removed `⋅¹` from line 598. Changed `\odot` to `\cdot\^2` in testing for double contraction (L487).
src/TensorValues/Operations.jl Outdated Show resolved Hide resolved
@fverdugo
Copy link
Member

We could probably consider also this notation...

¹,²,³,⁴,... # n-contractions
¹ == 
² == 

I guess we will need 4-contractions of 4-tensors, which are usually expressed as :: or 3-contractions, etc.

What do you think @fverdugo @Omega-xyZac ?

I dont know if these names are going to be interpreted as operators by the julia compiler

@zjwegert
Copy link
Contributor Author

We could probably consider also this notation...

¹,²,³,⁴,... # n-contractions
¹ == 
² == 

I guess we will need 4-contractions of 4-tensors, which are usually expressed as :: or 3-contractions, etc.
What do you think @fverdugo @Omega-xyZac ?

I dont know if these names are going to be interpreted as operators by the julia compiler

This seems to be working so far in my testing and other code.

@zjwegert zjwegert changed the title 3-tensor contractions with 2-tensor and 1-tensor 3-tensor contractions with 2-tensor and 1-tensor & 3-tensor zero method Oct 12, 2020
@zjwegert
Copy link
Contributor Author

#416 addressed here.

@zjwegert zjwegert changed the title 3-tensor contractions with 2-tensor and 1-tensor & 3-tensor zero method 3-Tensor Operations and Constructors Oct 12, 2020
Copy link
Member

@fverdugo fverdugo left a comment

Choose a reason for hiding this comment

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

I have added some more comments (see inline)

Note that ⋅² and are the same. I would only keep one version, perhaps ⋅². On the other hand, having double_contraction is fine since it is the plain asci version.

# a_ijpm = b_ijkl*c_klpm
@generated function double_contraction(a::A, b::B) where {A<:SymFourthOrderTensorValue{D},B<:SymFourthOrderTensorValue{D}} where D

Sym4TensorIndexing = [1111, 1121, 1131, 1122, 1132, 1133, 2111, 2121, 2131, 2122, 2132, 2133,
Copy link
Member

Choose a reason for hiding this comment

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

This is particular for 3D. Can you implement it in general?
If for performance reasons, you want to keep the 3D specialization then use the signature

function double_contraction(a::A, b::B) where {A<:SymFourthOrderTensorValue{3},B<:SymFourthOrderTensorValue{3}}

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Okay, perhaps ⋅². I will keep the 3D version since this should be fairly quick. I will have to have a think about how to do this generally.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This might do the trick:

@generated function double_contraction(a::SymFourthOrderTensorValue{D}, b::SymFourthOrderTensorValue{D}) where D
  str = ""
  for j in 1:D
    for i in j:D
      for m in 1:D
        for p in m:D
          s = ""
          for k in 1:D
            for l in 1:D
              s *= " a[$i,$j,$k,$l]*b[$k,$l,$p,$m] +"
            end
          end
          str *= s[1:(end-1)]*", "
        end
      end
    end
  end
  Meta.parse("SymFourthOrderTensorValue{D}($str)")
end

Need to test.

Copy link
Member

Choose a reason for hiding this comment

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

Yes, keep ⋅² so we use the same notation for all >1 contractions.

- Added generalised contraction.
- Removed ⊡.
@fverdugo fverdugo merged commit 9a09687 into gridap:master Oct 14, 2020
@fverdugo
Copy link
Member

@Omega-xyZac

can you include some comment on the new functionality in the https://github.com/gridap/Gridap.jl/blob/master/NEWS.md file ?

Please add the comments in a new section called [Unreleased], i.e.:

## [Unreleased]

## Added
- Your commend there...

Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants